Detecting subsystem symmetry protected topological order via entanglement entropy

TL;DR

This paper introduces a method to detect subsystem symmetry protected topological order through a constant correction to entanglement entropy, using tensor networks, and explores its robustness and extensions in various models.

Contribution

It provides an analytic and numerical framework to identify SSPT order via entanglement entropy corrections, including a new understanding of its persistence beyond the cluster phase.

Findings

The correction to the entanglement area law reliably detects SSPT order.

The SSPT order persists under certain external fields and interactions.

The correction behavior extends to three-dimensional cluster states.

Abstract

Subsystem symmetry protected topological (SSPT) order is a type of quantum order that is protected by symmetries acting on lower-dimensional subsystems of the entire system. In this paper, we show how SSPT order can be characterized and detected by a constant correction to the entanglement area law, similar to the topological entanglement entropy. Focusing on the paradigmatic two-dimensional cluster phase as an example, we use tensor network methods to give an analytic argument that almost all states in the phase exhibit the same correction to the area law, such that this correction may be used to reliably detect the SSPT order of the cluster phase. Based on this idea, we formulate a numerical method that uses tensor networks to extract this correction from ground-state wave functions. We use this method to study the fate of the SSPT order of the cluster state under various external fields and interactions, and find that the correction persists unless a phase transition is crossed, or the subsystem symmetry is explicitly broken. Surprisingly, these results uncover that the SSPT order of the cluster state persists beyond the cluster phase, thanks to a new type of subsystem time-reversal symmetry. Finally, we discuss the correction to the area law found in three-dimensional cluster states on different lattices, indicating rich behavior for general subsystem symmetries

Figures

Equations

$$S_A=a|\partial A|-\gamma +\dots$$

$$K_{x,y}|C⟩=|C⟩,$$

$${\mathcal{H}}_C=-\sum _{x,y}{K}_{x,y}.$$

$$U_v(c)$$

$$U_h(c)$$

$$S_A^{(\alpha )}\equiv S^{(\alpha )}({\rho }_A)=\frac{1}{1-\alpha }{\log }_2\mathrm{Tr}({\rho }_A^{\alpha }),$$

$$S_A=|A|-{\log }_2|G_A|.$$

$$H=\sqrt{\frac{1}{2}}\left(\begin{array}{cc}[r]1& 1\\ 1& -1\end{array}\right),$$

$$S_A^{(2)}(\theta )=N{\log }_{2}r-{\log }_{2}m+\dots ,$$

$$S_A^{(2)}(\pi )=\{\begin{array}{ll}2N-4& \text{mbox}if4\mid N\\ 2N-3& \text{mbox}if2\mid N\text{mbox}and4\nmid N\\ 2N-2& \text{mbox}if2\nmid N\end{array}.$$

$$A_2^i=B_2^i\otimes C_2^i.$$

$$\Pi =\frac{1+X^{\otimes N}}{2}$$

$$\mathbb{T}=\sum _{\mathbf{i}\in I}{\mathcal{A}}^{\mathbf{i}}\otimes {\mathcal{A}}^{\mathbf{i}\ast },$$

$$\mathcal{T}(\rho )=\sum _{\mathbf{i}\in I}{\mathcal{A}}^{\mathbf{i}}\rho {\mathcal{A}}^{\mathbf{i}†}$$

$$\mathcal{T}(\left[1\otimes X_k\right]\rho \left[1\otimes X_k\right])$$

$$=\sum _{\mathbf{i}\in I}{\mathcal{A}}^{\mathbf{i}}\left[1\otimes X_k\right]\rho \left[1\otimes X_k\right]{\mathcal{A}}^{\mathbf{i}†}$$

$$=\sum _{\mathbf{i}\in I_e}\left({\mathcal{B}}^{\mathbf{i}}\otimes {\mathcal{P}}^{\mathbf{i}}\Pi X_k\right)\rho \left({\mathcal{B}}^{\mathbf{i}}\otimes X_k\Pi {\mathcal{P}}^{\mathbf{i}}\right)$$

$$=\sum _{\mathbf{i}\in I_e}\left({\mathcal{B}}^{\mathbf{i}}\otimes X_k{\mathcal{P}}^{\mathbf{i}}\Pi \right)\rho \left({\mathcal{B}}^{\mathbf{i}}\otimes \Pi {\mathcal{P}}^{\mathbf{i}}{X}_k\right)$$

$$=\left[1\otimes X_k\right]\mathcal{T}(\rho )\left[1\otimes X_k\right].$$

$$\mathcal{T}(Z_k{Z}_{k+1}\rho Z_k{Z}_{k+1})=Z_k{Z}_{k+1}\mathcal{T}(\rho )Z_k{Z}_{k+1}.$$

$$\mathcal{T}(\rho )=X^{\otimes N}\mathcal{T}(\rho )=\mathcal{T}(\rho )X^{\otimes N}.$$

$$\left[{\sigma }_{R/L},X_k\right]$$

$$X^{\otimes N}{\sigma }_{R/L}$$

$$\mathcal{T}({\tilde{\sigma }}_R\otimes \Pi )$$

$$=\sum _{\mathbf{i}\in I}{\mathcal{B}}^{\mathbf{i}}{\tilde{\sigma }}_R{\mathcal{B}}^{\mathbf{i}†}\otimes {\mathcal{C}}^{\mathbf{i}}\Pi {\mathcal{C}}^{\mathbf{i}†}$$

$$=\sum _{\mathbf{i}\in I_e}{\mathcal{B}}^{\mathbf{i}}{\tilde{\sigma }}_R{\mathcal{B}}^{\mathbf{i}†}\otimes \Pi$$

$$\tilde{\mathcal{T}}({\tilde{\sigma }}_R):=\sum _{\mathbf{i}\in I_e}{\mathcal{B}}^{\mathbf{i}}{\tilde{\sigma }}_R{\mathcal{B}}^{\mathbf{i}†}={\tilde{\sigma }}_R.$$

$$S_A=S(\sigma )=S(\tilde{\sigma })+S(\Pi ),$$

$$S(\tilde{\sigma })\propto N$$

$$S_A^{(2)}=-{\log }_2\mathrm{Tr}({\sigma }_R{\sigma }_L^T{\sigma }_R{\sigma }_L^T)=-{\log }_2\mathrm{Tr}\left(Q^N\right),$$

$$S_A^{(2)}=N{\log }_2{\lambda }_0-\gamma +\dots ,$$

$$\mathcal{H}={\mathcal{H}}_C+\sum _{x,y}{h}_{x,y}^{\prime },$$

$$h_{x,y}^{\prime }=-h_X{X}_{x,y}.$$

$$h_{x,y}^{\prime }=-J_{ZZZZ}^{\diamond }{Z}_{x-1,y}{Z}_{x+1,y}{Z}_{x,y-1}{Z}_{x,y+1}.$$

Full text

The modern perspective of quantum phases of matter is based on global patterns of entanglement Chen et al. (2010); Wen (2017); Zeng et al. (2019). Two quantum states are said to lie in distinct (symmetry-protected) topological phases when they cannot be connected by (symmetric) local unitary evolution. As such, topological phases of matter can be characterized and detected by entanglement-based quantities. A well-known example is the topological entanglement entropy (TEE): For states with non-trivial topological order, the scaling of entanglement entropy exhibits a correction to the area law Hamma et al. (2005); Levin and Wen (2006); Kitaev and Preskill (2006); Eisert et al. (2010). That is, for a subregion $A$ of a lattice, the entropy for ground states of gapped, local Hamiltonians takes the general form,

[TABLE]

for some constants $a,\gamma$, where $\partial A$ is the boundary of region $A$ and the dots indicate terms that decay exponentially with $|A|$. The correction $\gamma$ takes a uniform non-zero value within non-trivial topological phases as a consequence of the global entanglement patterns. It can therefore be used to detect and characterize topological order and topological phase transitions analytically, numerically, and potentially even experimentally Papanikolaou et al. (2007); Hamma et al. (2008); Nussinov and Ortiz (2009); Isakov et al. (2011); Jiang et al. (2012); Flammia et al. (2009); Halász and Hamma (2012); Abanin and Demler (2012); Daley et al. (2012); Islam et al. (2015).

Recently, however, it has been observed that $\gamma$ may deviate from the expected value due to the presence of symmetry-protected topological (SPT) order localized around $\partial A$ Cano et al. (2015); Zou and Haah (2016); Santos et al. (2018); Devakul et al. (2018); Williamson et al. (2019); Schmitz et al. (2019). One setting in which this occurs is for states with subsystem SPT (SSPT) order Raussendorf et al. (2019); You et al. (2018a); Devakul et al. (2018, 2019). Such order is non-trivial only in the presence of subsystem symmetries, which are defined as symmetries that act on rigid lower-dimensional subsystems of the entire system. In cases where these symmetries act on 1D lines spanning a 2D lattice, one may find a non-zero value of $\gamma$ for regions $A$ whose boundaries lie parallel to these lines, despite the absence of topological order Zou and Haah (2016); Devakul et al. (2018); Williamson et al. (2019).

On one hand, these corrections may be viewed as obstructions to reliably extracting the TEE from ground states, and previous research has been focused on developing methods to get around this Zou and Haah (2016); Williamson et al. (2019). On the other hand, they suggest the possibility that SSPT order may also be characterized by corrections to the area law. This is an attractive prospect, as SSPT order has recently garnered notable interest in both the contexts of condensed matter, thanks to its relation to fracton order Vijay et al. (2016); Williamson (2016); Kubica and Yoshida (2018); You et al. (2018b); Shirley et al. (2019); Song et al. (2019); Ma et al. (2018), and quantum information, due to its use in measurement-based quantum computation Else et al. (2012a); Raussendorf et al. (2019); Stephen et al. (2019); Devakul and Williamson (2018). Thus, any tool that can be used to characterize and detect SSPT order would also have immediate impact in these areas.

Can $\gamma$ provide such a tool? Up to this point, SSPT order has been given as a sufficient condition for $\gamma\geq\log 2$ Devakul et al. (2018), but also an example of two states in the same SSPT phase with different non-zero values of $\gamma$ has been given Williamson et al. (2019). This suggests that SSPT phases of matter are not associated with a particular value of $\gamma$, and therefore that the precise value of $\gamma$ is not useful for the characterization of SSPT order. In this work, we demonstrate the contrary. We consider the 2D cluster phase, previously discussed in Refs. Raussendorf et al. (2019); Else et al. (2012a); You et al. (2018a); Devakul et al. (2018), and show that $\gamma$ takes the same value for all generic states in the phase, with deviations only occurring at fine-tuned points. This value of $\gamma$, which we refer to as the symmetry-protected entanglement entropy (SPEE), relates to the non-trivial symmetry fractionalization that occurs on the boundary of every state in the cluster phase. Therefore, the SPEE may be used to characterize SSPT phases of matter in the same way that TEE characterizes topological phases of matter.

We use this result as the basis for a new numerical technique to detect SSPT order in ground states of gapped local Hamiltonians. Namely, we show how one can straightforwardly extract $\gamma$ from a projected entangled pair state (PEPS) representation of a ground state, and then apply this method to various Hamiltonians obtained by perturbing the cluster state Hamiltonian. In accordance with our analytical arguments, we find that $\gamma$ can reliably and unambiguously detect SSPT order and SSPT phase transitions. Moreover, since $\gamma$ makes no reference to any particular symmetry of the system, it detects any and all SSPT order within a ground state, unlike usual SPT order parameters which must be defined with respect to a specific symmetry Bartlett et al. (2010); Haegeman et al. (2012); Pollmann and Turner (2012); Marvian (2017); You et al. (2018a). As a consequence of this, we discover a large region in which the SSPT order of the cluster state appears to persist despite the subsystem symmetries being explicitly broken. We make a preliminary analysis of this new phase of matter in terms of a new notion of subsystem time-reversal symmetry of the cluster state.

Together, our analytical and numerical results show that the SPEE is an effective tool to detect and characterize the SSPT order of the 2D cluster phase. Going beyond this, we also study 3D cluster states with different types of subsystem symmetries and calculate $\gamma$ in each case, observing distinct behaviours depending on the structure of the symmetries. We believe that our uniformity arguments for the 2D cluster phase will hold equally well for these 3D phases, and also other types of SSPT order. In an outlook, we discuss the implications of our results for measurement-based quantum computation, detection of fracton order, and the possible experimental observation of SSPT order.

This paper is organized as follows. In Sec. I, we define the cluster state and cluster phase, and show their relation to corrections to the area law. In Sec. II, we present an analytical argument that $\gamma$ is the same for almost every state in the cluster phase. Then, in Sec. III, we formulate our numerical method and use it to examine the SSPT order of the cluster state with various Hamiltonian fields and interactions added. In Sec. IV, we calculate the area law corrections for cluster states defined on various 3D lattices, and finally in Sec. V, we present our conclusions and future directions of work.

I Cluster phase and corrections to the area law

We will begin by reviewing the definition of the cluster state and the cluster phase. We will then calculate the entanglement entropy for a continuous one-parameter family of states in the cluster phase, and find that all states in this family have the same correction to the area law, except for a singular point where the correction is larger. We trace back this larger correction to extra symmetries of the state which do not generically hold in the cluster phase, thereby establishing the premise that all generic states in the cluster phase have the same correction to the area law.

Throughout this paper, we consider the geometry of a $45^{\circ}$ rotated 2D square lattice on an infinitely long cylinder with circumference $N$ (although similar results hold also for a finite cylinder or torus with dimensions much larger than the correlation length). The cylinder is bipartitioned into two semi-infinite subsystems, denoted $A$ and $B$. The 2D cluster state $|C\rangle$ Briegel and Raussendorf (2001); Raussendorf et al. (2003) can be defined on this lattice by placing a qubit in the $|+\rangle=\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)$ state on each lattice site, and applying the unitary $CZ=\mathrm{diag}(1,1,1,-1)$ between each pair of neighbouring qubits. It may equivalently be defined as the unique solution to the constraints,

[TABLE]

where $K_{x,y}=X_{x,y}Z_{x+1,y}Z_{x-1,y}Z_{x,y+1}Z_{x,y-1}$ and $X_{x,y}$ ($Z_{x,y},Y_{x,y}$) denotes the Pauli operator $\sigma^{X}$ ($\sigma^{Z},\sigma^{Y}$) acting on site $(x,y)$, see Fig. 1. This means that the cluster state is a stabilizer state. It is therefore the unique ground state of the following Hamiltonian:

[TABLE]

This Hamiltonian commutes with subsystem symmetries consisting of $X$ acting on every site along any diagonal line on the lattice (corresponding to vertical and horizontal lines along the cylinder):

[TABLE]

We set our origin such that $U_{v}(0)$ corresponds to the symmetry lying parallel to the boundary of $A$, see Fig. 1. The cluster state has SSPT order with respect to these symmetries Raussendorf et al. (2019); You et al. (2018a); Devakul et al. (2018), and we call the corresponding SSPT phase of matter the cluster phase Raussendorf et al. (2019).

This paper deals with calculating the entropy of entanglement of the reduced density matrix $\rho_{A}$, where the $A$ subsystem corresponds to the right half of the cylinder. For much of this paper we will make statements about the structure of $\rho_{A}$ directly, such that our claims hold for any $\alpha$-Rényi entropy

[TABLE]

including the Von Neumann entropy $S^{(1)}_{A}=-\mathrm{Tr}(\rho_{A}\log_{2}\rho_{A})$ obtained in the limit $\alpha\rightarrow 1$. However, we will sometimes focus on $S_{A}^{(2)}$ as it is most amicable to our numerical methods. It is also the most convenient to measure experimentally Abanin and Demler (2012); Daley et al. (2012); Islam et al. (2015). Whenever we make a statement that holds for all $\alpha$, we will simply use the notation $S$ or $S_{A}\equiv S(\rho_{A})$.

As the cluster state is a stabilizer state, $S_{A}$ may be straightforwardly calculated. Let $G$ be the group generated by all stabilizers $K_{x,y}$, and let $G_{A}\subset G$ be the group of all elements of $G$ which act non-trivially only on region $A$. Then we have the following equation Hamma et al. (2005):

[TABLE]

All stabilizers $K_{x,y}$ corresponding to lattice sites $(x,y)\in A-\partial A$ are contained in $G_{A}$. The product of all stabilizers along the boundary, which is precisely the line symmetry $U_{v}(0)$, is also in $G_{A}$, see Fig. 1. Hence we have $|G_{A}|=2^{|A|-|\partial A|+1}$, so $S_{A}=N-1$. We see that the SPEE takes the value $\gamma=1$ $(=\log_{2}2)$ for the cluster state, due to the subsystem symmetries forming non-local constraints lying along the boundary of $A$ Zou and Haah (2016); Williamson et al. (2019); Schmitz et al. (2019).

As a first venture away from the cluster state, we consider a family of states $|C(\theta)\rangle=U(\theta)|C\rangle$, where the circuit $U(\theta)$ is defined by acting with the two-body unitary $(H\otimes H)C\theta(H\otimes H)$ on every pair of neighbouring sites, with

[TABLE]

and $C\theta=\mathrm{diag}(1,1,1,e^{i\theta})$. $U(\theta)$ is diagonal in the local $X$-basis, hence it commutes with the subsystem symmetries of the cluster state, so $|C(\theta)\rangle$ is in the cluster phase for all $\theta$. We choose this family since it provides a smooth interpolation between the cluster state, $|C(0)\rangle$, and the state considered in Ref. Williamson et al. (2019), $|C(\pi)\rangle$, which was shown to display an enlarged value of $\gamma$.

For general $\theta$, $|C(\theta)\rangle$ is not a stabilizer state, so we need a different method to calculate its entropy. In order to calculate $S_{A}^{(2)}(\theta)$, we use the method of Ref. Zou and Haah (2016). Notice that $|C(\theta)\rangle$ may be created by a unitary circuit acting on a product state. Therefore, by applying unitaries on the $A$ and $B$ subsystems separately, which does not change the entropy, we may disentangle all qubits except those on a strip along the boundary. Thus the calculation of entropy for our 2D system is reduced to that of a 1D system with an extensive bipartitioning, and this may be easily computed using a transfer matrix method. Namely, we can construct a matrix $Q(\theta)$ that follows from the definition of $|C(\theta)\rangle$, such that $S_{A}^{(2)}(\theta)=-\log_{2}\mathrm{Tr}(Q(\theta)^{N})$, see Ref. Zou and Haah (2016) for more details. Let $\{\lambda_{k}\}_{k}$ be the set of eigenvalues of $Q(\theta)$ with maximum magnitude and write $\lambda_{k}=re^{i\phi_{k}}$. The entropy is then,

[TABLE]

where $m=\sum_{k}e^{iN\phi_{k}}$ and the terms contained in the dots decay exponentially in $N$. We can identify the constants $a$ and $\gamma$ from Eq. (1) with $\log_{2}r$ and $\log_{2}m$, respectively. Thus, the SPEE can be determined by examining the eigenvalues of $Q(\theta)$ with largest magnitude. If there are no non-positive eigenvalues with magnitude $r$, then the constant $m$ is simply the degeneracy of the largest eigenvalue. If such eigenvalues do exist, then $m$ can exhibit periodic $N$ dependence.

The spectrum of $Q(\theta)$ can be computed exactly 111This is possible because $Q(\theta)$ is simply a 256 $\times$ 256 matrix, which can be exactly diagonalized numerically, and the results are shown in Fig. 2. For $\theta\neq\pi$, there are two eigenvalues of largest magnitude, both of which are positive. Hence we have $m=2$ and $\gamma=1$. At $\theta=\pi$, there are 16 eigenvalues of largest magnitude. More precisely, we find that the non-zero eigenvalues of $Q(\pi)$ are $\frac{1}{4}$, $-\frac{1}{4}$, $\frac{i}{4}$, and $-\frac{i}{4}$ with degeneracies 8, 4, 2 and 2, respectively. In this case, we have $m=\left(8+4(-1)^{N}+2i^{N}+2(-i)^{N}\right)$, which gives,

[TABLE]

Therefore, the states $|C(\theta)\rangle$ have a SPEE $\gamma=1$ for all $\theta\neq\pi$, while $\gamma>1$ for $\theta=\pi$.

The enlarged SPEE for $\theta=\pi$ can be attributed to the fact that $|C(\pi)\rangle$ has many extra symmetries, aside from those in Eq. (4), which are not satisfied for $\theta\neq\pi$. Namely, $|C(\pi)\rangle$ is a stabilizer state satisfying $K^{\prime}_{x,y}|C(\pi)\rangle=|C(\pi)\rangle$ where the stabilizers $K^{\prime}_{x,y}$ are as pictured in Fig. 1. One can alternatively compute $S_{A}^{(2)}(\pi)$ using Eq. (6), and the results agree with Eq. (9) 222The enlarged SPEE comes from extra subsystem symmetries which appear in $G_{A}$. It would be interesting to investigate whether or not the enlarged SPEE would persist if these extra symmetries where preserved.. At the end of the next section, we will use the tensor network representation of $|C(\pi)\rangle$ to more clearly identify the mechanism through which these extra symmetries lead to a larger SPEE.

II Analytical argument for uniformity of the SPEE

The results of the previous section suggest that the SPEE $\gamma$ is uniform throughout the cluster phase, except for certain fine-tuned states with enhanced symmetries. Now, we will use tensor networks to give an analytical argument that all generic states in the cluster phase do indeed have the same SPEE. In fact, an elementary argument based purely on symmetry shows that the entanglement spectrum is $2^{N-1}$-fold degenerate, which implies that $S_{A}\geq N-1$, as described briefly in Appendix A. However, this does not tell us anything about the value of $\gamma$. The fact that $\gamma\geq 1$ was argued in Ref. Devakul et al. (2018) (using the result of Ref. Zou and Haah (2016)) under a similar set of assumptions as those used here 333Namely, they assumed that, throughout the phase, we can disentangle the degrees of freedom far from the cut using unitary operators localized on either side of it. This is essentially the same as assuming a finite bond-dimension PEPS representation, as we do here.. The purpose of the more careful argument of this section is to argue that $\gamma=1$ exactly. We begin by reviewing the tensor network characterization of the cluster phase given in Ref. Raussendorf et al. (2019), which forms the basis of our argument, and then we show how this can be used to constrain $\gamma$.

II.1 Tensor network description of the cluster phase

In Fig. 3, we define two PEPS representations of the cluster state. The first tensor, $C_{1}$, is the usual PEPS representation of the cluster state on the rotated square lattice. In Ref. Raussendorf et al. (2019), it was shown that every state in the cluster phase admits a PEPS representation on this lattice of the form,

[TABLE]

where $B_{1}$ is some tensor describing a projected entangled pair operator which commutes locally with $X$ as indicated. Throughout the cluster phase, $C_{1}$ stays the same, but $B_{1}$ varies 444We note that the construction in Ref. Raussendorf et al. (2019) gives PEPS tensors with exponentially large bond dimension, but this can be improved to give constant bond dimension when one takes advantage of the finite depth of the quantum circuits considered within.. In this section, we will find it more convenient to consider a different PEPS representation of the cluster state, denoted $C_{2}$, which has two qubits per unit cell and is defined on a non-rotated square lattice, see Fig. 3. Then the statement corresponding to Eq. (10) is,

[TABLE]

$B_{2}$ is again some tensor that varies through the phase. $X_{a/b}$ denotes $X$ acting on one of the two qubits in a unit cell, labelled $a/b$ as in Fig. 3.

If we contract the physical leg of any tensor $A$ with some state $|i\rangle\in\{|++\rangle,|+-\rangle,|-+\rangle,|--\rangle\}$, we obtain a tensor with 4 virtual legs, denoted $A^{i}$. In terms of these components, Eq. (11) reads,

[TABLE]

Thus, the virtual space of the tensor network decomposes into two subspaces, which are called the junk and protected subspaces, corresponding to the operators $B_{2}^{i}$ and $C_{2}^{i}$, respectively Else et al. (2012b). This decomposition is central to our arguments. It shows that the global pattern of entanglement that characterizes the cluster phase is that of the cluster state. An analogous statement is known to be true in certain 1D SPT phases, including that which contains the 1D cluster state Else et al. (2012b).

We now map our system onto a quasi-1D system by combining $2N$ spins around the cylinder into one larger block as shown in Fig. 1. The associated tensor $\mathcal{A}$ is obtained by combining $N$ tensors $A_{2}$ in a ring around the cylinder, see Fig. 4. If we let $\mathbf{i}=(i_{1},i_{2},\dots,i_{N})$ be an element of the index set $I$ labelling all $2^{2N}$ states in a block, we can again define the tensor components $\mathcal{A}^{\mathbf{i}}$. Due to the subsystem symmetries $U_{v}(c)$, the wavefunction of our state consists only of terms with an even number of $|-\rangle$ states in both columns in each block. Denoting this even-parity subset of states $I_{e}\subset I$ and the corresponding odd-parity subspace as $I_{o}\subset I$, we can therefore modify the tensor $\mathcal{A}$ by setting $\mathcal{A}^{\mathbf{i}}=0$, $\forall\mathbf{i}\in I_{o}$, without changing the state described by $\mathcal{A}$. This modification simplifies our description later on.

For the remaining $\mathbf{i}\in I_{e}$, Eq. (12) implies the decomposition $\mathcal{A}^{\mathbf{i}}=\mathcal{B}^{\mathbf{i}}\otimes\mathcal{C}^{\mathbf{i}}$, where $\mathcal{B}$ and $\mathcal{C}$ denote the blocked tensors living in the junk and protected subspaces, respectively. We can straightforwardly determine $\mathcal{C}^{\mathbf{i}}$ using the definition of $C_{2}$. We find that $\mathcal{C}^{\mathbf{i}}=\mathcal{P}^{\mathbf{i}}\Pi$, $\forall\mathbf{i}\in I_{e}$, where each $\mathcal{P}^{\mathbf{i}}$ is some tensor product of Pauli operators and

[TABLE]

is a rank-$2^{N-1}$ projection matrix with $X^{\otimes N}\equiv\prod_{k=1}^{N}X_{k}$. Therein, $X_{k}$ denotes $X$ acting on the $k$-th component of the $N$-component protected subspace. Importantly, the operators $\mathcal{P}^{\mathbf{i}}$ also satisfy $\mathcal{P}^{\mathbf{i}}\Pi=\Pi\mathcal{P}^{\mathbf{i}}$, for all $\mathbf{i}\in I_{e}$.

In Appendix A, we show how the non-trivial SSPT order of the cluster phase can be understood in terms of symmetry fractionalization on the virtual boundary of the PEPS. This follows directly from Eq. (11), and is the PEPS equivalent of the picture developed in Ref. Devakul et al. (2018).

II.2 Constraining $\gamma$ in the cluster phase

We will now use the above characterization of the cluster phase to show that the entanglement entropy $S_{A}$ can be decomposed into two parts, $S_{A}=S_{B}+S_{C}$, where $S_{C}$ is the entropy of the cluster state, and $S_{B}$ is essentially the entropy of the PEPS defined by the tensor $B^{i}_{2}$. Since we know $S_{C}=N-1$, and we will argue $S_{B}$ satisfies an area law with no correction for generic $B^{i}_{2}$, we will find that $S_{A}$ has a correction of $\gamma=1$ for generic states in the cluster phase.

First, we review the standard tensor network technique for calculating entanglement entropy Cirac et al. (2011). We define the transfer matrix

[TABLE]

where the star denotes complex conjugation. We can normalize our PEPS such that the largest eigenvalue of $\mathbb{T}$ is 1, and denote by $|R\rangle$ and $\langle L|$ the left and right fixed-points of $\mathbb{T}$, such that $\mathbb{T}|R\rangle=|R\rangle$ and $\langle L|\mathbb{T}=\langle L|$. For the following proof, we will find it more convenient to work with the associated quantum channel

[TABLE]

which is related to $\mathbb{T}$ by redefining the input and output legs of the tensor. We will refer to both $\mathbb{T}$ and $\mathcal{T}$ as the transfer matrix. The fixed points $|R\rangle$ and $\langle L|$ of $\mathbb{T}$ correspond to fixed points $\sigma_{R}$ and $\sigma_{L}$ of $\mathcal{T}$ and its adjoint $\mathcal{T}^{\dagger}=\sum_{\mathbf{i}}\mathcal{A}^{\mathbf{i}\dagger}\rho\mathcal{A}^{\mathbf{i}}$, respectively. See Fig. 4 for a graphical representation of these objects.

An important property of the cluster phase is that the fixed point of the transfer matrix is unique, with all other eigenvalues having magnitude less than unity. Specifically, a degenerate fixed point space of the transfer operator would manifest itself as long-range correlations along the cylinder axis – either for local operators, amounting to conventional long-range order Rispler et al. (2015), or for loop operators acting around the cylinder, amounting to topological order Schuch et al. (2013). Since any such order is absent in the cluster state (whose transfer operator has a unique fixed point), and thus by definition in the whole cluster phase, the fixed point of the transfer matrix must be unique in the whole cluster phase.

The reduced density matrix $\rho_{A}$ corresponding to half of the cylinder can be determined via a bulk-boundary correspondence for PEPS Cirac et al. (2011). Namely, we can write $\rho_{A}=V\sigma V^{\dagger}$ where $V$ is an isometry and $\sigma=\sqrt{\sigma_{L}^{T}}\sigma_{R}\sqrt{\sigma_{L}^{T}}$, with $T$ denoting the transpose. Thus, $S_{A}=S(\sigma)$. The rest of this section is dedicated to constraining $\sigma$ within the cluster phase.

To begin, note that, due to Eq. (12), $\mathcal{T}$ satisfies many symmetries throughout the entire cluster phase. First, we have,

[TABLE]

Therein, $X_{k}$ acts on the protected subspace, while $\mathbb{1}$ acts on the junk subspace. In the second equality, we substituted $\mathcal{A}^{\mathbf{i}}=\mathcal{B}^{\mathbf{i}}\otimes\mathcal{P}^{\mathbf{i}}\Pi$. In the third equality, we used the facts that $X_{k}\mathcal{P}^{\mathbf{i}}=\pm\mathcal{P}^{\mathbf{i}}X_{k}$ and $X_{k}\Pi=\Pi X_{k}$. In the same way, we can derive,

[TABLE]

Above, and in the rest of this section, we have omitted the identity operators acting on the junk subspace for notational simplicity. Finally, because $X^{\otimes N}\Pi=\Pi$, we also have the symmetries,

[TABLE]

All of the above symmetries hold also for $\mathcal{T}^{\dagger}$. These symmetries are a reflection of the symmetry fractionalization on the boundary that characterizes the cluster phase, as described in Appendix A.

Since the fixed-points are unique, they inherit these symmetries. For example, if $\sigma_{R}$ is the fixed-point of $\mathcal{T}$, Eq. (16) implies that $X_{k}\sigma_{R}X_{k}$ is as well. Since the fixed-point is unique up to a multiplicative constant, we have $X_{k}\sigma_{R}X_{k}\propto\sigma_{R}$. Since $\sigma_{R}$ is positive, and conjugation by $X_{k}$ preserves the trace, we must have $X_{k}\sigma_{R}X_{k}=\sigma_{R}$. Using similar arguments, we get all of the following symmetries:

[TABLE]

These symmetries completely constrain $\sigma_{R/L}$ on the protected subspace, with $\Pi$ being the unique solution to the constraints. Thus, the fixed points decompose across the junk and protected subspaces as $\sigma_{R/L}=\widetilde{\sigma}_{R/L}\otimes\Pi$ for some unknown states $\widetilde{\sigma}_{R/L}$. This gives $\sigma=\widetilde{\sigma}\otimes\Pi$, where $\widetilde{\sigma}=\sqrt{\widetilde{\sigma}_{L}^{T}}\widetilde{\sigma}_{R}\sqrt{\widetilde{\sigma}_{L}^{T}}$.

Having constrained the fixed points in the protected subspace, what can we now say about the unconstrained part $\widetilde{\sigma}$? It turns out that $\widetilde{\sigma}_{R/L}$ are themselves fixed points of a certain transfer operator. Indeed, we have,

[TABLE]

where we used the facts that $\mathcal{C}^{\mathbf{i}}=0$, $\forall\mathbf{i}\in I_{o}$, and $[\mathcal{C}^{\mathbf{i}},\Pi]=0$, $\forall\mathbf{i}\in I_{e}$. Since $\mathcal{T}(\widetilde{\sigma}_{R}\otimes\Pi)=\widetilde{\sigma}_{R}\otimes\Pi$ by definition, this gives

[TABLE]

Furthermore, Eq. (20) implies that every eigenvector $\widetilde{\rho}$ of $\widetilde{\mathcal{T}}$ yields an eigenvector $\widetilde{\rho}\otimes\Pi$ of $\mathcal{T}$ with the same eigenvalue. Therefore, $\widetilde{\sigma}_{R}$ is the unique fixed-point of $\widetilde{\mathcal{T}}$, since the fixed-point of $\mathcal{T}$ is unique. Similarly, $\widetilde{\sigma}_{L}$ is the unique fixed-point of $\widetilde{\mathcal{T}}^{\dagger}$.

Since $\sigma=\widetilde{\sigma}\otimes\Pi$, the entanglement entropy decomposes into two components,

[TABLE]

where $S(\Pi)=N-1$ is the entropy of the cluster state. What remains is to understand the entropy contribution $S(\widetilde{\sigma})$. We will argue that $S(\widetilde{\sigma})$ generically satisfies an area law with no correction. To this end, consider a generic tensor $A_{2}^{i}=B_{2}^{i}\otimes C_{2}^{i}$ which is in the cluster phase. Therefore, it can be connected to the cluster state via a smooth path $A_{2}^{i}(\theta)=B_{2}^{i}(\theta)\otimes C_{2}^{i}$ without closing the gap of the transfer operator $\mathcal{T}$. This implies that $B_{2}^{i}(\theta)$ smoothly connects $B_{2}^{i}$ to a trivial tensor, where the corresponding transfer operator $\widetilde{\mathcal{T}}$ must remain gapped along the path as well. Up to the restriction $\mathbf{i}\in I_{e}$, $\widetilde{\mathcal{T}}$ is thus the transfer operator of a system in the trivial phase. However, we expect that such a global parity constraint will only affect the entropy if the system either has topological order – which we have ruled out – or physical symmetries in the basis of the constraint, which would require fine-tuning of the $B_{2}^{i}$. Thus, we expect the fixed point of $\widetilde{\mathcal{T}}$ to have the structure of a generic fixed point in the trivial phase, which does not exhibit long-range correlations and thus no constant correction to the area law.

Overall, this reasoning implies that, for generic points in the cluster phase, the entropy $S(\widetilde{\sigma})$ should scale as

[TABLE]

with corrections vanishing as $N\to\infty$, whereas constant corrections are expected only at fine-tuned points with additional symmetries. This is confirmed by numerical study of generic tensors $B_{2}^{i}$ up to bond dimension $4$. We thus find that, for a generic point in the cluster phase, the entanglement entropy should scale like $S_{A}=aN-1$ for some constant $a$.

In light of these arguments, let us now reconsider why $|C(\pi)\rangle$ from Sec. I has a correction $\gamma>1$. For this state, one can write a PEPS in the form of Eq. (11), where $B^{i}$ is the same tensor as $C^{i}$, up to applying $H$ on each physical leg. Therefore, in addition to Eqs. (16-18), $\mathcal{T}$ has just as many extra symmetries which act non-trivially on the junk subspace. These symmetries serve to constrain $\widetilde{\sigma}_{R/L}$, leading to the increased value of $\gamma$. Such symmetries are generically not present in the cluster phase, and indeed disappear for all $\theta\neq\pi$, reflected by the generic correction $\gamma=1$.

III Numerical detection of SSPT order

Above, we have argued that $\gamma=1$ within the entire cluster phase, except at certain fine-tuned points of enhanced symmetry. In this section, we perform numerical calculations of $\gamma$ in ground states. The motivation of this is twofold. First, we would like to substantiate our analytical arguments. By considering known models, we will confirm that $\gamma=1$ within the cluster phase, and $\gamma=0$ in the trivial phase. Second, once we are confident that $\gamma=1$ indicates the presence of SSPT order, we can use it as a probe for phase transitions. We will add various fields and interactions onto the cluster Hamiltonian, with phase transitions into the trivial phase indicated by $\gamma$ suddenly dropping from 1 to 0. This study will lead us to the discovery of a new phase of matter beyond the cluster phase.

III.1 Description of the algorithm

We use the following method to calculate $\gamma$ in ground states. Given a Hamiltonian $H$, we use infinite projected entangled pair state (iPEPS) techniques to obtain a PEPS approximation to the ground state Vanderstraeten et al. (2016). We then use this PEPS to construct the transfer matrix as in Eq. (15), and use an infinite matrix product state (iMPS) algorithm to obtain matrix product operator (MPO) approximations for the fixed-points $\sigma_{R}$ and $\sigma_{L}$. Using this, we can write

[TABLE]

where $Q$ is a matrix obtained from the MPO tensors describing $\sigma_{R}$ and $\sigma_{L}^{T}$, see Fig. 4 555We remark that the matrix $Q$ is not the same object as the matrix $Q(\theta)$ introduced in Sec. I, although they are closely related. For all states encountered in this section, the leading eigenvalue of $Q$ will either be unique or two-fold degenerate. Denote the two largest eigenvalues of $Q$ (which may or may not be degenerate) as $\lambda_{0}$ and $\lambda_{1}$. In analogy with Eq. (8), we have

[TABLE]

where $\gamma=1$ if $\lambda_{0}=\lambda_{1}$, and $\gamma=0$ otherwise. Therefore, we can determine $\gamma$ from the ratio $\lambda_{1}/\lambda_{0}$. We note that this ratio is closely related to the replica correlation length that was introduced in Ref. Zou and Haah (2016) to study spurious corrections to the area law. A similar method to compute entanglement entropy in PEPS was also proposed in Ref. Poilblanc et al. (2016).

III.2 Stability of the cluster phase

Now we use the above method to investigate the fate of the SPEE under various local fields or interactions added to the cluster Hamiltonian, one at a time. The goal is to support our argument that the value $\gamma=1$ characterizes the cluster phase, meaning that this value persists unless the subsystem symmetries are explicitly broken or a phase transition is crossed.

We consider a Hamiltonian of the form:

[TABLE]

for several different choices of $h^{\prime}_{x,y}$, each with different symmetries. In Fig. 5, we plot $\lambda_{1}/\lambda_{0}$ as a function of the strength of the perturbation $h^{\prime}_{x,y}$. The important findings are the following:

(i) When a small symmetry respecting perturbation is added to the Hamiltonian, the SPEE keeps the value $\gamma=1$.

(ii) When the strength of the symmetry respecting perturbations are increased, we encounter phase transitions indicated by a sudden drop to $\gamma=0$. The location of the phase transitions agree with known results where available.

(iii) Adding a term which does not commute with the subsystem symmetries removes the SPEE for any finite coupling strength.

Overall, these results support our analytical results and confirm that the SPEE is a useful probe for detecting SSPT order and phase transitions. We will now discuss each choice of $h^{\prime}_{x,y}$ in more detail.

III.2.1 Symmetry respecting terms

We first consider adding terms which commute with the subsystem symmetries. For such terms, we expect that the value of the SPEE will remain constant up until an SSPT phase transition is reached, at which point it should should jump to 0. The first term we consider is a simple local $X$-field:

[TABLE]

The model of a cluster Hamiltonian with added $X$-field been studied previously in Refs. Doherty and Bartlett (2009); Kalis et al. (2012); Orús et al. (2013); You et al. (2018a). For $h_{X}\rightarrow\infty$, the model becomes a trivial paramagnet without SSPT order, so there must be a phase transition at some value of $h_{X}$. The model is self dual under the unitary circuit $U_{CZ}:=U(\pi)$ (see Sec. I), so this phase transition should occur at $h_{X}=1$. Indeed, via a non-local duality transformation, this model can be mapped to the so-called Xu-Moore model with a transverse field, which is known to have a first order phase transition at $h_{X}=1$ Xu and Moore (2004); Doherty and Bartlett (2009); You et al. (2018a). In agreement with these facts, we find that $\lambda_{1}/\lambda_{0}$ changes from 1 discontinuously at the transition point, such that $\gamma$ jumps to 0 at this point, thereby correctly detecting the phase transition.

The next term we consider is a four-body interaction,

[TABLE]

This is the minimal term that contains $Z$ operators yet still commutes with all subsystem symmetries. Interestingly, this model behaves in the opposite way to first one, in that it is mapped to the Xu-Moore model under $U_{CZ}$, while it is self-dual under the same non-local duality transformation. These facts predict a phase transition at $J^{\diamond}_{ZZZZ}=1$, and this is again in agreement with the behaviour of $\lambda_{1}/\lambda_{0}$.

The final symmetry-preserving term we consider is a nearest-neighbour Ising interaction,

[TABLE]

To the best of our knowledge, this model has yet to be studied in the literature. Via the same non-local duality transformation mentioned earlier, it may be mapped onto two Xu-Moore models coupled by Ising interactions. The SPEE disappears at $J_{XX}=0.5$, suggesting that there is an SSPT phase transition into the symmetry-breaking phase at this point.

We also verified that simultaneously adding all three symmetry respecting terms with random small couplings does not change the SPEE. This shows that the behaviour observed above is not a consequence of the specific Hamiltonians considered, and is instead generic behaviour.

III.2.2 Symmetry breaking terms

We now turn our attention to terms which anti-commute with the subsystem symmetries, and therefore explicitly break them. For these terms, the SSPT order should be destroyed for any finite coupling strength, so $\lambda_{1}/\lambda_{0}$ should decrease from 1 immediately. The simplest symmetry-breaking term is a $Z$-field,

[TABLE]

As predicted, we find $\lambda_{1}/\lambda_{0}<1$ for any value of $h_{Z}$. Furthermore, this model is mapped to a trivial Hamiltonian under the action of $U_{CZ}$. Therefore, there is no phase transition for any value of $h_{Z}$, and this is consistent with the smooth decay of $\lambda_{1}/\lambda_{0}$ that we observe.

The next term we consider is a next-nearest neighbour interaction,

[TABLE]

While this term anti-commutes with the subsystem symmetries, it commutes with the global “checkerboard” $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ symmetry of the cluster state, which is generated by applying $X$ to all even or all odd spins on the square lattice. The cluster state has “weak” 2D SPT order under this global symmetry group, which is defined by a 2D SPT order that is only non-trivial in the presence of translational invariance Chen et al. (2013); Miller and Miyake (2016). Therefore, it is important to confirm that the non-zero value of the SPEE is due to the subsystem symmetries, and not the global symmetries alone. Indeed, we find $\lambda_{1}/\lambda_{0}<1$ for any finite value of $J^{\times}_{ZZ}$, indicating that the global symmetries are not sufficient to protect the SPEE. As opposed to the previous case, we observe some singular behaviour of $\lambda_{1}/\lambda_{0}$ as $J^{\times}_{ZZ}$ is increased. This is explained by the fact that, under $U_{CZ}$, this model is mapped to two decoupled 2D quantum Ising models, which undergo a second order phase transition into a symmetry-breaking phase at $J^{\times}_{ZZ}\approx 0.3285$. Thus, the behaviour of $\lambda_{1}/\lambda_{0}$ also helps to detect non-SSPT phase transitions.

III.3 Beyond the cluster phase — Time reversal symmetry

We now consider adding the term,

[TABLE]

This term anti-commutes with the subsystem symmetries, as with Eq. (30). Therefore, according to our findings in the previous subsection, we would expect any finite value of $h_{Y}$ to remove the SPEE. This turns out not to be the case, as seen in Fig. 6, and this section is devoted to understanding this behaviour.

It was shown in Ref. Stephen et al. (2019) that the cluster state also has non-trivial SSPT order protected by fractal subsystem symmetries which are composed of tensor products of $Y$ operators. Since the model is trivial when $h_{Y}\rightarrow\infty$, there should be a phase transition into the trivial phase for some $h_{Y}$. Since the model is self-dual under $U_{CZ}$ followed by applying $S=\mathrm{diag}(1,i)$ on every site, this transition should occur at $h_{Y}=1$. Along the $y$-axes of Fig. 6, we see good evidence that this is the case. For example, the magnetization along the $Y$ axis jumps discontinuously at this point.

What is surprising is that this transition is also detected by the SPEE, which remains equal to 1 up until $h_{Y}=1$, even though the $Y$-field does not commute with the subsystem symmetries. It is tempting to attribute the SPEE to the fractal symmetries of the cluster state, which do commute with the $Y$-field, but this is not the case. Consider simultaneously adding $X$ and $Y$ magnetic fields to the cluster state,

[TABLE]

We see in Fig. 6 that the SPEE still persists in a large region. This cannot be due to the $X$ linelike symmetries of Eq. (4), nor the $Y$ fractal symmetries, since both are explicitly broken by the added fields. Therefore, it appears that we have discovered a new SSPT phase which contains at least part of cluster phase, and which is also accompanied by a non-zero SPEE.

A similar phenomenon occurs in 1D. The 1D cluster state has SPT order under a global $\mathbb{Z}_{2}\times\mathbb{Z}_{2}$ symmetry, generated by applying $X$ to all even or all odd spins on the chain. Similar to the current scenario, this SPT order is stable under adding a $Y$-field, despite the fact that this is a symmetry-breaking term. In this case, the resolution is that the 1D cluster state also has time-reversal symmetry Verresen et al. (2017). This symmetry also protects the SPT order of the cluster state, and commutes with the $Y$-field.

Time-reversal turns out to be the solution here as well, although it takes an unusual form. In Ref. You et al. (2018a), the authors defined subsystem time-reversal symmetries as subsystem unitary symmetries, such as the linelike symmetries considered here, followed by global time-reversal. But it is easy to see that this also does not work, since the global time-reversal flips the sign of all $Y$-fields, not just those lying along a given line of subsystem symmetry. We therefore need a different notion of subsystem time-reversal, in which we enact time-reversal locally only on those sites of the lattice on which the subsystem symmetry acts non-trivially.

Ref. Chen and Vishwanath (2015) suggested one way to implement time-reversal locally on a tensor network. Globally, time-reversal symmetry acts on a quantum state by the action of a unitary operator (here, the subsystem symmetries) combined with complex conjugation of the wavefunction. With tensor networks, the wavefunction is divided into local tensors. Thus, one can define local complex conjugation at a given site by conjugating only the tensor associated to that site. More precisely, we can define operators $\mathcal{K}_{x,y}$ which act on the PEPS by conjugating the tensor at site $(x,y)$ only. Note that, to ensure that the proper notion of locality is used, we use the tensor network defined on a rotated square lattice with a single-qubit unit cell, shown in Fig. 3.

With this, we define our subsystem time-reversal as,

[TABLE]

Is this the correct symmetry to describe the phase of matter observed in Fig. 6? To answer this, we consider a state described by a tensor of the following form,

[TABLE]

where $C_{1}$ is the cluster tensor, and $B_{T}^{*}$ is the complex conjugate of $B_{T}$. The state described by a tensor of this form is not generally symmetric under $U_{h,v}(c)$, but it is symmetric under $U^{T}_{h,v}(c)$ (note that the cluster state tensor $C_{1}$ is real, and hence unaffected by $\mathcal{K}$). Furthermore, if we pick a random bond dimension 4 tensor $A_{T}$ of the form given in Eq. (35), then we can confirm numerically that the resulting PEPS does indeed have $\gamma=1$.

We conjecture that ground states in the non-trivial phase in Fig. 6 can be captured by tensors of the form Eq. (35), in the same way that the cluster phase is captured by tensors of the form Eq. (10). In Appendix B, we argue that this is indeed the case by using the framework of perturbation theory in PEPS Vanderstraeten et al. (2017). For example, to first order in perturbation theory, we can approximate the ground states of the model described by Eqs. (26,33) with the following PEPS tensor,

[TABLE]

where the matrix $\Lambda$ is defined as,

[TABLE]

Notice that $\Lambda$ satisfies $X\Lambda X=\Lambda^{*}$, such that $A$ is of the form given in Eq. (35). Furthermore, one can confirm that the fixed-point $\sigma$ of the associated transfer matrix is exactly $\Pi$, independent of the values of $h_{X}$ of $h_{Y}$, so we have $\gamma=1$.

A similar argument holds to all orders in perturbation theory. Therefore, at least up to the phase transition where perturbation theory breaks down, the ground states in the non-trivial phase of Fig. 6 can be captured by PEPS of the form Eq. (35), which we believe will generically have $\gamma=1$, using similar arguments as in Sec. II. On the other hand, if we also add a $Z$ field to the Hamiltonian, the tensor obtained from perturbation theory no longer has the form of Eq. (35) and, accordingly, the SPEE disappears.

Overall, our understanding of the phase of matter in Fig. 6 comes from three steps: (1) The subsystem time reversal symmetries $U^{T}_{h,v}(c)$ are the relevant symmetries that protect the SSPT order and (2) states in the same phase as the cluster state under $U^{T}_{h,v}(c)$ can be captured by PEPS of the form Eq. (35), which (3) generically have a SPEE of $\gamma=1$. We leave a more rigorous confirmation of steps (1)-(3) to future work.

III.4 Remarks on the numerical method

We finish this section by commenting on some relevant features of our numerical method. To begin, we briefly compare our method of computing $\gamma$ to the more common approach, which involves calculating ground states and their respective entropies $S_{A}$ for several system sizes $N$, and then making a linear fit to extract $\gamma$ Jiang et al. (2012). This method is relatively costly, in that it requires the determination of several ground states, and it is also sensitive to finite size effects. On the other hand, our method can obtain $\gamma$ from one PEPS tensor, and works directly in the thermodynamic limit, such that finite size effects are minimized.

Second, we see from Figs. 5 and 6 that the SPEE can serve as a very good probe for SSPT phase transitions, with the phase boundary marked by a clear discontinuity. In some cases, the phase transitions observed here may also be detected by discontinuities in local magnetizations. For example, the phase transition induced by an added $X$-field can be detected by $\langle m_{X}\rangle$, as seen in Fig. 6. However, Fig. 6 also shows that neither $\langle m_{X}\rangle$ nor $\langle m_{Y}\rangle$ can resolve the entire phase diagram alone. This is to be expected, as SPT phases cannot be completely detected by local order parameters. The SPEE, on the other hand, is inherently non-local, and can resolve the entire phase diagram. Nevertheless, the SPEE is still comparably simple to calculate given a PEPS ground state, so it is a genuinely useful tool to detect SSPT phase transitions.

Finally, we note a distinct property that the SPEE has as a tool to detect SSPT order. Usual quantities used to detect SPT order, such as string order parameters, are explicitly defined in terms of a certain symmetry Bartlett et al. (2010); Haegeman et al. (2012); Pollmann and Turner (2012); Marvian (2017); You et al. (2018a). Such quantities therefore can only determine whether a state has SPT order with respect to this symmetry. In contrast, the SPEE is “symmetry-agnostic”, meaning that it is not defined with respect to any particular symmetry of the state 666The geometry of the bipartition we use to calculate entropy is informed by the geometry of the subsystem symmetries, but this is the only way in which the symmetry enters the definition of the SPEE. This is the reason that we were able to detect the phase of matter in Fig. 6 without a prior understanding of the relevant symmetry. The fact that a symmetry-agnostic quantity like the SPEE can characterize SSPT order reflects its fundamental differences from standard SPT order Schmitz et al. (2019).

IV 3D cluster states

We have now seen that the 2D cluster phase can be detected and characterized by a value $\gamma=1$ of the SPEE. In this section, we move to investigating the SPEE in 3D systems. Going to 3D allows us to consider different types of subsystem symmetry, and also brings our discussion closer to fracton order, which exists only in dimension 3 and higher. We will consider cluster states defined on different 3D lattices, each with a different notion of subsystem symmetry, and each leading to a different type of area law correction. In light of our results for the 2D cluster phase, we expect the observed behaviour to hold throughout the corresponding SSPT phases.

We start by reviewing how to define cluster states on lattices other than the 2D square lattice, also known as graph states. Consider a graph $\mathcal{G}=(V,E)$ consisting of vertices $v\in V$ and edges $e\in E$. To define the corresponding graph state $|\mathcal{G}\rangle$, we place a qubit in the $|+\rangle$ state on every vertex, and then act with $CZ$ on every pair of vertices connected by an edge,

[TABLE]

$|\mathcal{G}\rangle$ is the ground state the following Hamiltonian:

[TABLE]

where $N(v)\subset V$ is the neighbourhood of $v$, and contains all vertices $u$ where are connected to $v$ by an edge.

To compute entanglement entropy, we need only consider an effective 2D system lying along the boundary of the partition, as discussed in Sec. I. For these 2D boundary systems, we can use Eq. (6) to calculate $S_{A}$. In what follows, we will consider 3D cubic lattices with dimensions $L\times L\times L$ (such that $|\partial A|=L^{2}$) and periodic boundary conditions in all directions. Similar statements can be made for open boundary conditions after the addition of boundary terms in the Hamiltonians. We will consider a biparition whose boundary is perpendicular to a coordinate axis, and all symmetries discussed will be tensor products of Pauli-$X$ operators. We note that graph state symmetries made of Pauli-$Y$ operators can also exist, but with drastically different, e.g. fractal, geometry as shown in Ref. Stephen et al. (2019).

Our results are summarized in Fig. 7. We will now briefly discuss each model.

(i) bcc lattice. This model has subsystem symmetries in the form of lines moving in the $x,y,z$ directions of the lattice. It was previously studied in Ref. You et al. (2018a), where it was also shown that it has non-trivial SSPT order under these linelike symmetries. The 2D system living on the boundary of the bipartition is exactly the 2D cluster state, where the $A$ and $B$ subsystems correspond to the $A$ and $B$ sublattices. Therefore, the $2L-1$ linelike symmetries of the cluster state generate $G_{A}$ in Eq. (6), and we get a lower dimensional “perimeter law” correction to the area law.

(ii) fcc lattice. The second model has linelike symmetries on the face qubits, but only planar symmetries on the vertex qubits. The 2D boundary system has only a single global symmetry in $G_{A}$, so there is a constant correction rather than a perimeter law correction. Thus, the planar symmetries seem to dictate the physics of this model, rather than the linelike symmetries. This example, along with the previous, suggest that SSPT phases in dimension $D$ with $k$-dimensional subsystem symmetries are associated with a correction to the area law that scales like $L^{D-k-1}$.

(iii) RBH lattice. The third model, first introduced by Raussendorf, Bravyi, and Harrington (RBH) Raussendorf (2005), was originally used in the context of fault-tolerant quantum computation Raussendorf et al. (2006, 2007). It was later shown in Ref. Roberts et al. (2017) to have SPT order protected by so-called “1-form” symmetries Yoshida (2016). Such symmetries are distinct from subsystem symmetries because they are deformable. While subsystem symmetries have rigid geometry like lines or planes, 1-form symmetries in a 3D system can live on any closed 2D surface. This follows from local symmetries generated by applying $X$ to every face qubit in a single unit cell (or every edge qubit connected to a given vertex). In particular, this model also has planar symmetries as in the fcc lattice. Despite this, the boundary state is trivial, consisting of entangled pairs between the two subsystems, so there is no correction to the area law. The reason for this is as follows: For the previous two cases, as well as the 2D cluster state, the correction emerges because there is a non-local stabilizer that lives in one subsystem and is composed of local stabilizers that have support in both subsystems. For the case of 1-form symmetries, this is not the case: the non-local stabilizers can be decomposed into local stabilizers that are also contained in one subsystem. This behaviour was also noticed in a 2D example (which is actually the same state emerging on the boundary of the fcc lattice), in Ref. Williamson et al. (2019). Therefore, 1-form symmetries are likely not associated with area law corrections in general.

V Outlook

In this paper, we have shown that the 2D cluster phase can be characterized by a uniform correction $\gamma=1$ to the area law which we have called the symmetry-protected entanglement entropy (SPEE), and we used this to construct a new numerical technique to detect SSPT order. The SPEE is relatively easy to calculate given a ground state PEPS tensor, hence it is an effective tool to detect SSPT order and SSPT phase transitions. A next step would be to check if the results presented in this paper are unique to the 2D cluster phase, or are generally true for all non-trivial SSPT phases. For example, Ref. Devakul et al. (2018) showed that the SSPT phases protected by the linelike symmetries considered here can be sorted into 8 different equivalence classes, and constructed representative states for each class. It is also easy to generalize the cluster state to different symmetry groups as in Ref. You et al. (2018a) and Ref. Brell (2015), for example. Also, as shown here, there is rich behaviour of the correction for 3D cluster states. In all of these cases, it would be interesting to see whether there is a uniform value of the SPEE within the corresponding SSPT phases.

Using our method, we uncovered the surprising result that the SSPT order of the cluster state is preserved under the addition of local magnetic fields pointing anywhere in the $X$-$Y$ plane, as indicated by the SPEE. This indicates that the cluster phase defined previously Raussendorf et al. (2019); You et al. (2018a); Devakul et al. (2018) is part of a larger, more robust SSPT phase of matter. We attributed this larger phase of matter to subsystem time-reversal symmetry, although an in-depth understanding is still missing. In particular, the physical meaning of our notion of subsystem time-reversal is unclear. Interestingly, the newly discovered phase (‘T’) contains at least part of both the cluster phase (‘X’) and the fractal SSPT phase (‘Y’) containing the cluster state Stephen et al. (2019). This suggests that the phases ‘X’ and ‘Y’, whose symmetries differ drastically in geometry, can potentially be unified by considering the anti-unitary symmetries of the phase ‘T’.

This newly discovered phase may also have implications in terms of quantum computation. It is known that phases ‘X’ and ‘Y’ are computationally universal Raussendorf et al. (2019); Stephen et al. (2019), meaning that every state within these phases may be used as a resource for universal measurement-based quantum computation. If the phase ‘T’ is also computationally universal, this would mean that the computational power of the cluster state is robust to a much larger class of perturbations than previously thought. In general, it would be worthwhile to understand whether there is a clear link between the SPEE and the usefulness of a state as a computational resource, and if so, how well the SPEE compares to other figures of merit used to determine the computational usefulness of, e.g., perturbed cluster states Doherty and Bartlett (2009); Kalis et al. (2012); Orús et al. (2013).

In 3D, certain models with subsystem symmetries are dual to models with fracton order Vijay et al. (2016); Williamson (2016); Kubica and Yoshida (2018); You et al. (2018b); Shirley et al. (2019); Song et al. (2019). This duality is realized by gauging the subsystem symmetries. In particular, if this gauging procedure is applied to a state with SSPT order, the resulting fracton model can be twisted You et al. (2018b), in the same way that gauging a 2D SPT leads to a twisted Toric Code, i.e. Double Semion model Levin and Gu (2012). It is therefore plausible to think that a transition between different fracton orders could be dual to a transition between different SSPT orders, which could in turn be detected by the SPEE, or a 3D generalization thereof.

Finally, we briefly comment on the possibility to measure the SPEE experimentally. Measuring the topological entanglement entropy would be a direct way of observing the presence of topological order in a many-body state. However, this is experimentally daunting; a principle barrier being the difficulty of creating topological states in the first place. On the other hand, cluster states are relatively easy to make in optical lattices with nearest-neighbour interactions Mandel et al. (2003). Furthermore, the second Rényi entropy can be measured in optical lattices by interfering two identical copies of a ground state and performing local parity measurements Abanin and Demler (2012); Daley et al. (2012), and this has been done to observe area laws in 1D systems Islam et al. (2015). Thus, measuring the SPEE of the cluster phase seems feasible with current or near-term technologies. This would serve as a first route to verifying non-trivial quantum order via entanglement entropy.

Acknowledgements.

DTS was supported by a fellowship from the Natural Sciences and Engineering Research Council of Canada (NSERC). This work has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme through the ERC Starting Grant WASCOSYS (No. 636201), and by the Deutsche Forschungsgemeinschaft (DFG) under Germany’s Excellence Strategy (EXC-2111 – 390814868).

Appendix A Subsystem symmetry fractionalization in PEPS

In this appendix, we show how the non-trivial SSPT order of the cluster phase may be seen analytically by observing symmetry fractionalization on the PEPS boundary. Symmetry fractionalization, referring to the notion that the symmetry of an SPT ordered state acts in a fundamentally different way on the boundary of the system compared to in the bulk, is the basic principle of SPT order. For example, the boundary of 1D SPTs transforms protectively under the symmetry group Pollmann et al. (2010); Chen et al. (2011a); Schuch et al. (2011), while the symmetry of 2D SPTs acts non-locally on the boundary Chen et al. (2011b). Understanding how the symmetry acts on the boundary is the key to characterizing and classifying SPT order Chen et al. (2013); Else and Nayak (2014). The case is the same for SSPT order, as discussed in Refs. You et al. (2018a); Devakul et al. (2018). For example, it was observed that, in the cluster phase, two parallel lines of subsystem symmetry which commute in the bulk may not commute on the boundary. This observation led to a classification of SSPT phases Devakul et al. (2018).

The framework of PEPS gives a natural way to understand the boundary of a many-body state. Specifically, for a PEPS defined on a region with boundaries, the uncontracted virtual legs at the boundary of the region allow to parametrize the Hilbert space for both the physical boundary excitations Yang et al. (2014) and the entanglement spectrum of a bipartition Cirac et al. (2011) on equal footing, and allow for a direct understanding of the way in which physical symmetries act on the boundary and entanglement.

In Fig. 8, we show how the subsystem symmetries act on these boundary degrees of freedom throughout the entire cluster phase. The relations pictured therein can be straightforwardly verified for the cluster state by using the definition of the tensor $C_{2}$. Looking just at the top edge of Fig. 8, we see that the boundary action of a given symmetry line anti-commutes with those of the neighbouring lines to the left and right (note that the $b$ qubit on each site is to the right of the $a$ qubit, cf. Fig. 3), and commutes with all others. This projective representation has irreducible dimension $2^{N-1}$. This same projective symmetry representation appears on the boundary of every state in the cluster phase, where the Pauli operators in the virtual space are acting on the protected subspace and the action on the junk subspace is trivial. One way to see this is to use Eq. (11). A more elementary argument starts by adopting a quasi-1D picture of our 2D system by wrapping it onto a long cylinder. Then, the fact that the projective representation on the boundary cannot change without crossing a phase transition follows from the standard arguments for 1D SPT phases Pollmann et al. (2010); Chen et al. (2011a); Schuch et al. (2011).

This symmetry fractionalization is a direct indication of the non-trivial SSPT order of the cluster phase. It can be used to derive, for example, the extensive ground state degeneracy of SSPT ordered systems with open boundaries You et al. (2018a), as well as a $2^{N-1}$-fold degeneracy in the entanglement spectrum on a cylinder of circumference $N$ Pollmann et al. (2010). The latter implies a bound on the entropy, $S_{A}\geq N-1$. Furthermore, our proof in Sec. II is based on these symmetries, which appear as symmetries of the transfer matrix in Eqs. (16)-(17). Therefore, the value $\gamma=1$ of the SPEE is a result of the way in which the subsystem symmetries fractionalize on the boundary.

Fig. 8 is the PEPS equivalent of the picture developed in Ref. Devakul et al. (2018), where the action of subsystem symmetries on a physical boundary are considered, as opposed to the virtual boundary considered here. The advantage of PEPS is that we can clearly see the same anti-commutation for all states in the cluster phase, and we are not restricted to models with zero correlation length. Different SSPT phases under the same subsystem symmetries would correspond to different anti-commutation relations on the virtual boundary. It should therefore be possible to reproduce the complete classification of Ref. Devakul et al. (2018) in the PEPS picture.

Appendix B Perturbation theory in PEPS

In this appendix, we use perturbation theory in PEPS to show that states in the non-trivial phase of Fig. 6 may be well approximated by tensors of the form Eq. (35). Consider a small perturbation away from the cluster state

[TABLE]

where $h^{\prime}_{x,y}=-h_{X}X_{x,y}-h_{Y}Y_{x,y}$ with $h_{X},h_{Y}\ll 1$. Throughout this appendix, we will label sites by a single index $j=(x,y)$ to condense notation. Defining

[TABLE]

the ground state, to first order in perturbation, is given by

[TABLE]

for

[TABLE]

In going from Eq. (43) to (44), we have used that (adopting the notation of Sec. IV),

[TABLE]

and similarly $\braket{C}{Y_{j}}{C}=0\quad\forall j$, since

[TABLE]

To go from Eq. (44) to (45), we exploited the fact that

[TABLE]

which implies that for all $j$, $X_{j}\ket{C}$ is an eigenstate of $(\mathcal{H}_{C}-E_{C})^{-1}$ with inverse energy 1/8 (note that $j\notin N(j))$. The same statement holds for $Y_{j}\ket{C}$ with inverse energy 1/10. The first order perturbation is implemented by acting on the cluster state by acting with $\Lambda$ on every site, cf. Eq. (36). Note that $\Lambda$ fulfills $X\Lambda X=\Lambda^{*}$, by virtue of only consisting of Pauli-$X$ and $Y$ matrices with real coefficients.

To reach higher orders in perturbation theory, we use the exponential perturbation theory developed in Ref. Vanderstraeten et al. (2017), which requires us to evaluate $\exp{(V)}|C\rangle$. Therefore, we must evaluate terms like $V^{2}\ket{C},V^{3}\ket{C}$, and so on. We can use the machinery of Vanderstraeten et al. (2017) to write down PEPOs enacting the perturbations with increasing bond dimension for any given order. The extra bond dimension is needed to account for the fact that higher powers of $V$ will introduce products of $X$s and $Y$s with different prefactors, depending on the geometry of the operator. For example, the second order contribution to the ground state is given by

[TABLE]

where the first four terms sum over nearest neighbours, next-nearest neighbours, next-next-nearest neighbours, and all other pairings, respectively. Together with the first order correction, one (not necessarily optimal) PEPO that implements this correction is given by a local tensor with bond dimension $D=11$ and non-zero entries (the order of the virtual legs is north, west, south, east) given by

[TABLE]

[TABLE]

which fullfils $XB_{ijkl}X=B^{*}_{ijkl}$ for all $i,j,k,l$. The tensor $B$ therefore commutes locally with $X$, up to complex conjugation, as in Eq. (35).

We will now argue why this holds to all orders in perturbation theory: As long as the action of $V^{n}\ket{C}$ on the cluster state can be expressed as sums of Pauli strings of $X$s and $Y$s with real coefficients, one can always find a PEPO $B$ for $\exp(V)$ which fulfills $XB_{ijkl}X=B^{*}_{ijkl}$. In other words, we want to show that

[TABLE]

for all $n$, where $P_{\mathbf{a}\mathbf{b}}=\bigotimes_{k}X^{a_{k}}Y^{b_{k}}$, $\mathbf{a}$ and $\mathbf{b}$ are bitstrings on the lattice and the $c_{\mathbf{a}\mathbf{b}}$ are real coefficients. We have seen explicitly that Eq. (59) holds for $n=1,2$. Now, assuming Eq. (59) holds for some $n$, we have,

[TABLE]

where $c^{\prime}_{\mathbf{a}\mathbf{b}}$ are also real. Therein, we have used that the perturbation $h^{\prime}$ consists of Pauli-$X$s and $Y$s with real coefficients and the fact that each product of Pauli matrices either maps the cluster state to an exact excitation or to itself, making $(\mathcal{H}_{C}-E_{C})^{-1}(\mathbb{1}-\ket{C}\bra{C})$ act simply as a multiplication of each of the $c_{\mathbf{a}\mathbf{b}}$ by a real number. So Eq. (59) holds for all $n$ by induction.

Therefore, it is clear that, at any given order, the corresponding PEPO-tensor $B$ describing the perturbation will again act as a real linear combination of $\mathbb{1},X$, $Y$ and $XY$ for any given virtual state, hence the resulting PEPS tensor has the form of Eq. (35). We conclude that there must exist a region in $(h_{X},h_{Y})$ space around the cluster point in which the ground state is accurately described by a PEPO fulfilling Eq. (35) acting on the cluster state.