Entanglement Spectra of Stabilizer Codes: A Window into Gapped Quantum Phases of Matter

TL;DR

This paper investigates the entanglement spectra of stabilizer codes to understand gapped quantum phases, revealing universal non-local features and the influence of subsystem symmetries, especially in fracton and SSPT phases.

Contribution

It develops a framework for computing entanglement features of gapped topological phases in stabilizer models, including fracton and SSPT phases, highlighting the role of non-local surface stabilizers.

Findings

Non-local surface stabilizers encode universal entanglement features.

Subsystem symmetries induce geometric dependence in entanglement spectra.

Edge-entanglement correspondence extends to 3D fracton models.

Abstract

The entanglement spectrum (ES) provides a barometer of quantum entanglement and encodes physical information beyond that contained in the entanglement entropy. In this paper, we explore the ES of stabilizer codes, which furnish exactly solvable models for a plethora of gapped quantum phases of matter. Studying the ES for stabilizer Hamiltonians in the presence of arbitrary weak local perturbations thus allows us to develop a general framework within which the entanglement features of gapped topological phases can be computed and contrasted. In particular, we study models harboring fracton order, both type-I and type-II, and compare the resulting ES with that of both conventional topological order and of (strong) subsystem symmetry protected topological (SSPT) states. We find that non-local surface stabilizers (NLSS), a set of symmetries of the Hamiltonian which form on the boundary of the entanglement cut, act as purveyors of universal non-local features appearing in the entanglement spectrum. While in conventional topological orders and fracton orders, the NLSS retain a form of topological invariance with respect to the entanglement cut, subsystem symmetric systems---fracton and SSPT phases---additionally show a non-trivial geometric dependence on the entanglement cut, corresponding to the subsystem symmetry. This sheds further light on the interplay between geometric and topological effects in fracton phases of matter and demonstrates that strong SSPT phases harbour a measure of quasi-local entanglement beyond that encountered in conventional SPT phases. We further show that a version of the edge-entanglement correspondence, established earlier for gapped two-dimensional topological phases, also holds for gapped three-dimensional fracton models.

Figures

a of Ref. Schmitz ([2018](#bib.bib81)) has two independent versions given by exchanging the configurations in the three $[111]$ layers of the triangular lattice over which the pattern is periodic; see the reference for more details. For example, two are implied by the topological constraints given by the product over all stabilizers of a single type, which enforces

Equations

$$O_s|\psi ⟩=|\psi ⟩,$$

$$H_{\mathcal{S}}=-\sum _s{J}_s{O}_s,(J_s>0)$$

$$g(\mathbf{n})=\prod _{i=1}^{d_G}{O}_i^{n_i}.$$

$${\mathcal{P}}^{\mathbf{k}}=\frac{1}{|G|}\sum _{g(\mathbf{n})\in G}(-1{)}^{\mathbf{k}\cdot \mathbf{n}}g(\mathbf{n}),$$

$$g(\mathbf{n}){\mathcal{P}}^{\mathbf{k}}=(-1{)}^{\mathbf{k}\cdot \mathbf{n}}{\mathcal{P}}^{\mathbf{k}}.$$

$$H=H_{\mathcal{S}}+\lambda V,$$

$$V=\sum _q\sum _{i=1}^3{\xi }_q^i{X}_q^i$$

$$H_{\text{CSS}}=-\sum O_s^X-\sum O_s^Z,$$

$${\mathcal{P}}^{\mathbf{k}}=|\mathbf{k}⟩⟨\mathbf{k}|,$$

$$E_{\mathbf{k}}=2\Vert \mathbf{k}\Vert -d_G,$$

$${\rho }_A={\text{Tr}}_B|\mathbf{k}⟩⟨\mathbf{k}|=2^{-s_A}{\mathcal{P}}_A^{{\mathbf{k}}_A},$$

$$s_A=|A|-d_{G_A}.$$

$$s_B=|B|-d_{G_B}-d_l=s_A.$$

$$\mu =\min (d_{\partial }-s_A-s_B),$$

$$G_{\text{NLSS}}=G_{\partial }\cap (G_A\oplus G_B),$$

$$g_{\partial }=g_A{g}_B,$$

$$\prod _{i\in F_A}{O}_i^A=(g_{\partial }{)}_A.$$

$$g_{\partial }=g_B,$$

$$|\mathbf{k}⟩=2^{-\frac{s}{2}}\sum _{{\mathbf{l}}_{\partial }}(-1{)}^{{\mathbf{l}}_{\partial }\mathbf{Mk}}{|{\mathbf{k}}_A,{\mathbf{l}}_{\partial }⟩}_A\otimes {|{\mathbf{k}}_B,{\mathbf{l}}_{\partial }\oplus {\mathbf{k}}_{\partial }⟩}_B,$$

$$(O_i{)}_A{|{\mathbf{k}}_A,{\mathbf{l}}_{\partial }⟩}_A=(-1{)}^{({\mathbf{l}}_{\partial }{)}_i}{|{\mathbf{k}}_A,{\mathbf{l}}_{\partial }⟩}_A,$$

$$\mathcal{L}=-\sum _{\mathbf{k},\mathbf{l}}\frac{V_{\mathbf{kl}}[E_{\mathbf{k}}\ne E_{\mathbf{l}}]}{E_{\mathbf{k}}-E_{\mathbf{l}}+0}|\mathbf{k}⟩⟨\mathbf{l}|,$$

$$H_{\text{ent}}\simeq -{\left(\frac{\partial {\tilde{\rho }}_A}{\partial \lambda }\right)}_{\lambda =0},$$

$$P^{\mathbf{k}}{X}_q^i=X_q^i{P}^{\mathbf{k}\oplus {\mathbf{p}}_q^i},$$

$$|⟨\mathbf{k}⟩X_q^i\mathbf{l}{|}^2=⟨\mathbf{k}⟩X_q^i{P}^{\mathbf{l}}{X}_q^i\mathbf{k}={\delta }_{\mathbf{k}\oplus \mathbf{l},{\mathbf{p}}_q^i}.$$

$$⟨\mathbf{k}⟩X_q^i\mathbf{l}=(-1{)}^{\mathbf{k}\oplus {\tilde{\mathbf{p}}}_q^i}{\delta }_{\mathbf{k}\oplus \mathbf{l},{\mathbf{p}}_q^i}.$$

$$V_{\mathbf{kl}}=⟨\mathbf{k}⟩V\mathbf{l}=\sum _q\sum _{i=1}^3(-1{)}^{\mathbf{k}\oplus {\tilde{\mathbf{p}}}_q^i}{\xi }_q^i{\delta }_{\mathbf{k}\oplus \mathbf{l},{\mathbf{p}}_q^i}.$$

$$E_{\mathbf{k}}-E_{\mathbf{k}\oplus {\mathbf{p}}_q^i}=2(\Vert \mathbf{k}\oplus {\mathbf{p}}_q^i\Vert -\Vert \mathbf{k}\Vert )=2(\overline{\mathbf{k}}\cdot {\mathbf{p}}_q^i-\mathbf{k}\cdot {\mathbf{p}}_q^i),$$

$$\mathcal{L}=\sum _{\mathbf{k}}\sum _q\sum _{i=1}^3{C}_{\mathbf{k}}^{q_i}|\mathbf{k}⟩⟨\mathbf{k}\oplus {\mathbf{p}}_q^i|,$$

$$C_{\mathbf{k}}^{q_i}=\frac{1}{2}\frac{(-1{)}^{\mathbf{k}\cdot {\tilde{\mathbf{p}}}_q^i}{\xi }_q^i[\mathbf{k}\cdot {\mathbf{p}}_q^i\ne \bar{\mathbf{k}}\cdot {\mathbf{p}}_q^i]}{\overline{\mathbf{k}}\cdot {\mathbf{p}}_q^i-\mathbf{k}\cdot {\mathbf{p}}_q^i+0}.$$

$$C_{\mathbf{k}\oplus {\mathbf{p}}_r^j}^{q_i}=C_{\mathbf{k}}^{q_i},$$

Full text

1. I Introduction

2. II Entanglement Spectra of Stabilizer Code Hamiltonians

3. II.1 The Stabilizer Formalism

4. II.2 Eigenstates and Entanglement Spectrum of Unperturbed Stabilizer Codes

5. II.3 Perturbed Density Matrix and Entanglement Hamiltonian

6. II.4 Entanglement Hamiltonian for different topological sectors

7. III Entanglement Spectra for Long-Range Entangled States

8. III.1 $\mathbb{Z}_{2}$ Topological Order

9. III.1.1 $d=2$ Toric Code/Wen-plaquette Model

10. III.1.2 $d\geq 3$ Toric Code

11. III.2 Type-I Fracton Order: X-cube model

12. III.3 Type-II Fracton Order: Haah’s Cubic Code

13. III.4 Edge-ES Correspondence

14. IV Quasi-Local Entanglement in SSPT Phases

15. V Discussion: Role of NLSS

16. VI Conclusions and Outlook

17. A Unitary Perturbation Theory

18. B Canonical Basis for the Pauli Group

19. C Recoverable Information for the $[111]$ cut in Haah’s code

I Introduction

The study of interacting many-body quantum phases of matter has been greatly influenced and aided by developments in the field of quantum information; indeed, it is now widely appreciated that quantum entanglement plays a vital role in the characterisation of a wide variety of zero temperature complex many-body systems. Entanglement has proven a particularly efficacious tool in the study of gapped quantum phases, especially in lower spatial dimensions ($d<3$) where it has provided a classification of gapped 1d phases Chen et al. (2010); Fidkowski and Kitaev (2011) and a powerful diagnostic for topological order in $d=2$ spatial dimensions Kitaev and Preskill (2006); Levin and Wen (2006). Gapped topologically ordered phases of matter are believed to be entirely characterised by the universal properties of their ground states, and a defining feature of these phases, which evade description in terms of any local order parameter, is the pattern of long range entanglement (LRE) in their ground state(s) (see e.g. Ref. Grover et al. (2013) for a review). A potent and oft utilised probe for the presence of LRE is the entanglement entropy—for topologically ordered systems in two spatial dimensions, the entanglement entropy contains a universal subleading ‘constant’ term, known as the “topological entanglement entropy (TEE),” which partially characterises such phases and is intimately linked to their topological quantum field theory (TQFT) description.

The entanglement spectrum (ES), introduced by Li and Haldane in the context of fractional quantum Hall fluids Li and Haldane (2008), provides a more general probe of quantum entanglement and encodes physical information beyond that contained in the entanglement entropy. For a given bipartition of a system into two regions $A$ and $B$, the “entanglement Hamiltonian” $H_{A}$ of its ground state $\ket{\Psi}$ is defined through $\rho_{A}\equiv e^{-H_{A}}/Z$, where $Z=\text{Tr}\,e^{-H_{A}}$ and $\rho_{A}=\text{Tr}_{B}\ket{\Psi}\!\bra{\Psi}$ is the reduced density matrix defined on region $A$. Since the entanglement cut mimics a physical boundary for the system, $H_{A}$ can be crudely understood as the Hamiltonian for a physical edge of the system; based on this observation, Li and Haldane conjectured that the ES contains universal information about the low-energy boundary excitations i.e., the ground state wave function in the bulk encodes information about dynamics at the edge.

Following the original proposal, the entanglement spectrum has been widely used to identify and distinguish quantum phases of matter, and has been especially fruitful in characterising gapped topological phases. While much of the initial work focused on chiral topological orders in 2d, including fractional quantum Hall fluids Regnault et al. (2009); Papić et al. (2011); Chandran et al. (2011); Dubail et al. (2012); Cano et al. (2015); Qi et al. (2012), this was later extended to include symmetry protected topological (SPT) states Swingle and Senthil (2012); Prodan et al. (2010); Pollmann et al. (2010); Turner et al. (2010); Fidkowski (2010); Alba et al. (2012); Choo et al. (2018) as well as topologically ordered states with gapped boundaries Ho et al. (2015, 2017); Koch-Janusz et al. (2017); Luo et al. (2018). In all cases, a key result is the existence of a correspondence between the low-lying spectrum of the physical edge states and the low-lying entanglement spectrum, often referred to as an edge-ES correspondence.

In three spatial dimensions, fracton order (see Ref. Nandkishore and Hermele (2018) for a review) has emerged as a new platform for realising long-range entanglement in gapped quantum phases111Gapless fracton order remains an equally active area of research, where symmetric tensor gauge theories have emerged as a powerful formalism within which to encapsulate much of the fracton phenomena Pretko (2017a, b); Prem et al. (2018); Gromov (2017); Ma et al. (2018a); Bulmash and Barkeshli (2018)—we will not discuss these here., exhibiting features both familiar and distinct from those encountered in conventional topological order. Familiar features include the presence of long-range entanglement, locally indistinguishable ground states on non-trivial manifolds, and topologically charged excitations which cannot be created locally. The crucial distinguishing feature of these phases—discovered first in a series of exactly solvable models Chamon (2005); Haah (2011); Yoshida (2013); Vijay et al. (2015, 2016)—which has engendered much recent activity Williamson (2016); Slagle and Kim (2017a); Prem et al. (2017); Hsieh and Halász (2017); Pretko and Radzihovsky (2018); Pai and Pretko (2018); Shirley et al. (2018); Song et al. (2018); Pai et al. (2018); Pretko and Radzihovsky (2018); Kumar and Potter (2018); You and von Oppen (2018), is that the mobility of certain excitations (fractons) is strictly verboten, while that of other excitations (sub-dimensional particles) is restricted along sub-dimensional manifolds of the three dimensional lattice. Unlike topologically ordered phases, whose ground state degeneracy depends only on the topology of the underlying manifold, for fracton phases the number of ground states grows sub-extensively, revealing their sensitivity to the underlying geometry. In this precise sense, fracton phases differ from topologically ordered phases in that they do not admit a low-energy description as a TQFT; indeed, a series of recent works have emphasised the geometric nature of fracton order Slagle and Kim (2017b, 2018); Shirley et al. (2018); Prem et al. (2018); Slagle et al. (2018a); Yan (2018); Slagle et al. (2018b); Gromov (2018).

Concurrently, the field of subsystem symmetry protected topological (SSPT) phases You et al. (2018); Devakul et al. (2019, 2018a); Kubica and Yoshida (2018); Stephen et al. (2018); Williamson et al. (2018a) has emerged as a close relative of fracton order, with a large number of fracton models obtained through gauging the subsystem symmetry of various SSPT phases Vijay et al. (2016); You et al. (2018, 2018); Shirley et al. (2018a); Williamson et al. (2018b). Subsystem symmetries, alternatively referred to as “gauge-like” symmetries Nussinov and Ortiz (2009), act on rigid subsystems which cannot be deformed, such as rigid planes of the cubic lattice. Besides their deep connection to fracton phases, it has also been realised that $d=2$ SSPT phases protected by rigid line-like symmetries can act as a universal resource for measurement-based quantum computation Else et al. (2012); Devakul and Williamson (2018); Stephen et al. (2018).

Given the theoretical discovery of these novel gapped phases, it is natural to ask whether measures of entanglement, such as the entanglement entropy or entanglement spectrum, exhibit features distinct from those seen in SPT phases or those with topological order. Early work in this direction has mostly focused on understanding the EE of these phases, with their ES remaining largely uncharted Shi and Lu (2018); He et al. (2018); Ma et al. (2018b); Schmitz et al. (2018); Shirley et al. (2018b); Williamson et al. (2018a); Yan (2018). In this paper, we aim to fill this lacuna by studying the entanglement spectra of stabilizer codes, which encode a large class of gapped quantum phases of matter. As such, stabilizer codes provide a complementary language to that of TQFTs, furnishing zero correlation length Hamiltonians which are sums of commuting projectors and which describe an exactly solvable point within some phase. Importantly, stabilizer codes provide a universal language within which we can compare and contrast the entanglement structure of fracton and SSPT phases with that of topologically ordered states.

Specifically, in this paper, we develop a general procedure for deriving the ES for the ground state of a stabilizer code Hamiltonian in the presence of arbitrary weak, local perturbations. Typically, when calculating the TEE in the ground state of some stabilizer Hamiltonian, one first derives the EE for the fixed point Hamiltonian and then argues that any contributions to the TEE are invariant under local perturbations. However, perturbations play a crucial role in revealing the structure of the ES and must be taken into account from the beginning. Indeed, while the ES is flat for any unperturbed stabilizer code, it presents universal features in the presence of weak, local perturbations for all models considered here.

As examples of topological order, we consider the $d=2,3$ toric code models and find that their entanglement Hamiltonians (EH) universally map onto $\mathbb{Z}_{2}$ invariant $d=1,2$ Ising models acting on effective spin-$1/2$ variables. As a consequence of this, we recover a version of the edge-ES correspondence first established in Ref. Ho et al. (2015). We then consider the X-cube model Vijay et al. (2016) and Haah’s cubic code Haah (2011) as examples of type-I and type-II fracton orders respectively and find that the geometric nature of these phases is manifest in their ES as well. That is, the EH for these models can be mapped onto an effective subsystem symmetric Ising-like model only for those entanglement cuts consistent with the planar (fractal) subsystem symmetries of the X-cube (cubic code). Thus, we show that the ES serves as a clear entanglement measure distinguishing fracton order from topological order, adding to the existing diagnostics for fracton order Devakul et al. (2018b); Prem et al. (2018); Schmitz et al. (2018). We also provide strong evidence for a correspondence between the low-lying ES and the low-lying spectrum of physical edge states, extending the validity of the edge-ES correspondence to gapped fracton phases. Finally, we consider the $d=2$ cluster state as our example of a strong SSPT phase and argue that the so-called “spurious” contributions to the TEE, found in Ref. Williamson et al. (2018a), are in fact evidence of quasi-local entanglement in such phases. In other words, we find that the entanglement structure of this system harbours features redolent of LRE systems which, despite not being topological in nature, remain robust against all subsystem symmetric perturbations. Thus, we establish SSPT phases as lying between conventional SRE SPT phases and LRE topological orders, thereby expanding the dictionary of possible patterns of entanglement in gapped phases.

The balance of this paper is organised as follows: in Section II we start by reviewing some basic notions of the stabilizer formalism pertinent to this paper. In Section II.2, we derive the flat entanglement spectra of any stabilizers’ eigenstate and discuss some consequences, including a review of recoverable information and emergent Gauss’ laws, as well as of the Schmidt decomposition for the stabilizers’ eigenstates. This is followed by the perturbative analysis in Section II.3, where we arrive at the central result of the paper: a general formula for the entanglement Hamiltonian of a stabilizer ground state in the presence of an arbitrary weak, local perturbation for any bipartition of the qubits $(A,B)$, with $B=A^{c}$. In Section III, we apply our general formalism to gapped LRE phases: the $d=2,3$ toric codes, the X-cube model, and Haah’s cubic code, with the latter two being examples of type-I and type-II fracton models, respectively. Section IV applies our general method to a strong-SSPT model, namely the $d=2$ cluster state, where we discuss signatures of subsystem symmetries in the entanglement and introduce the notion of quasi-local entanglement. Section V summarises the relation between all models considered herein and discusses a conjecture which extends our results to low-energy excited states. Lastly, we make our final remarks in Section VI.

II Entanglement Spectra of Stabilizer Code Hamiltonians

We begin this section by reviewing some basics regarding the stabilizer formalism before deriving general results for the structure of quantum entanglement in these systems. In certain cases, we rely upon results previously established in Refs. Haah (2013); Ma et al. (2018b); Schmitz et al. (2018), and provide details only when required for the remainder. After reviewing the requisite background, we present a general method for deriving the entanglement spectrum of the ground state of a stabilizer Hamiltonian in the presence of arbitrary perturbations, allowing us to compare and contrast a myriad of phases within the same framework. From the ES, one can also extract the entanglement entropy and recoverable information Schmitz et al. (2018); however, as the EE has been studied for both fracton order and SSPTs before, we only comment on this briefly.

II.1 The Stabilizer Formalism

The stabilizer formalism provides a common language within which a large class of quantum many-body systems and quantum error-correcting codes can be efficiently described (see Ref. Brown et al. (2016) for a review). It provides a unifying framework for studying several interacting quantum many-body systems, as certain exactly solvable points within some gapped phase can be described in terms of a Hamiltonian which is the sum of commuting Pauli operators—stabilizers. Each of the models considered here will admit such a description.

The models under consideration here consist of a set of $N$ local qubits (spin-1/2 degrees of freedom) living on the edges or vertices of a simple graph, with $q$ qubits per each edge or vertex, such that $N=q|V|$, where $V$ is the graph set. Here, $|\cdots|$ denotes the size of a finite set. The Hilbert space of the system $\mathcal{H}\simeq\mathbb{C}_{2}^{\otimes N}$ is given by the product space of all qubits.

The Pauli group $P$ acting on $N$ qubits is defined as the set of all Pauli operators—with individual qubit Pauli operators denoted by $\{X,Y,Z\}$—acting on the qubits, modulo any phase of $\pm 1,\pm i$. The stabilizer set $\mathcal{S}\subset P$ is a subset of the Pauli group comprised of mutually commuting operators, which satisfies $|\mathcal{S}|\geq N$ (there are at least as many stabilizers as qubits), $-I_{\mathcal{H}}\not\in\mathcal{S}$ (all of $\mathcal{S}$ can have a positive eigenvalue), and $\text{Supp}(\mathcal{S})=V$ (each element of $V$ is acted upon non-trivially by at least one stabilizer in $\mathcal{S}$). The set of states $\ket{\psi}$ which are stabilized by $\mathcal{S}=\{O_{s}\}$:

[TABLE]

for all $O_{s}\in\mathcal{S}$, form the ground state manifold for the stabilizer code Hamiltonian

[TABLE]

i.e., states invariant under the action of elements of $\mathcal{S}$ span the ground state subspace of the Hamiltonian $H_{\mathcal{S}}$.

Since all members of $\mathcal{S}$ commute, stabilizers in $\mathcal{S}$ multiplicatively generate an Abelian group $G=\{\prod_{s\in F}O_{s}:F\in\mathbb{P}[\mathcal{S}]\}$, where the power set of $\mathcal{S}$, denoted $\mathbb{P}[\mathcal{S}]$, is the set of all subsets of $\mathcal{S}$. Not all stabilizers are independent, so $\mathcal{S}$ may over-determine $G$. We refer to any $C\subseteq\mathcal{S}$ such that $\prod_{s\in C}O_{s}=I$ as a constraint. Following Haah (2013); Ma et al. (2018b); Schmitz et al. (2018), let $d_{G}=\log_{2}|G|$ and $\{O_{i}\}_{i\leq d_{G}}\subseteq\mathcal{S}$ be a complete, independent generating set for $G$. Elements $g\in G$ can hence be labelled by a binary vector $\mathbf{n}=(n_{1},n_{2},\dots,n_{d_{G}})\in\{0,1\}^{d_{G}}$ through

[TABLE]

With every vector $\mathbf{k}\in\{0,1\}^{d_{G}}$, we associate a projection operator

[TABLE]

where $\mathbf{k}\cdot\mathbf{n}$ is the binary dot product of $\mathbf{k}$ and $\mathbf{n}$. It is easy to check that $\left(\mathcal{P}^{\mathbf{k}}\right)^{2}=\mathcal{P}^{\mathbf{k}}$, all $\mathcal{P}^{\mathbf{k}}$ mutually commute, and that

[TABLE]

This implies that $\{\mathcal{P}^{\mathbf{k}}\}$ are the projection operators onto the simultaneous eigenstates of all of $G$ and, since $H_{\mathcal{S}}$ is a sum over a (over)complete generating set for $G$, onto the energy eigenvalue subspaces as labelled by $\mathbf{k}$. Clearly, $\mathbf{k}$ are the quantum numbers labeling excitations, where $\mathbf{k}=0$ labels the ground state manifold.

The projectors $\mathcal{P}^{\mathbf{k}}$ are pure state projectors if and only if $d_{G}=N$. In most cases under consideration, $d_{G}<N$ if the system is on a topologically non-trivial manifold, reflecting that all $\mathbf{k}$ eigenstate manifolds—including the ground state manifold—of the Hamiltonian (2) are degenerate on that manifold222This is a defining feature of systems with long-range entanglement; for SPTs, on the other hand, there is a unique ground state on arbitrary manifolds and, correspondingly, $d_{G}=N$ for those systems.. To see this, one can take the trace of $\mathcal{P}^{\mathbf{k}}$ and use the fact that tr$(g(\mathbf{n}))=N\delta_{\mathbf{n},0}$ to find that the dimension of the $\mathbf{k}$ eigenvalue manifold is $2^{d_{l}}$ for $d_{l}=N-d_{G}$. We refer to this type of degeneracy as a topological degeneracy. To obtain pure state projectors, we can complete $G$ by adding $d_{l}$ Pauli operators to our generating set, such that these operators are mutually commuting and also commute with all elements $g\in G$. These $d_{l}$ operators, which preserve the ground state manifold (also referred to as the code-space in the context of error correction) but act non-trivially on it, are referred to as the logical operators of the stabilizer code, a reflection of the fact that a stabilizer code encodes $d_{l}$ logical qubits in the ground space. Formally, the set of logical operators of a code is specified by $L\coloneqq\mathcal{C}(G)/G$, where $\mathcal{C}(G)$ is the centralizer of $G$.

In this paper, we study phases of matter captured by stabilizer code Hamiltonians of the form

[TABLE]

where $H_{\mathcal{S}}$ is given by Eq. (2) with $J_{s}=1$ for all $s$ and where

[TABLE]

describes a perturbation to $H_{\mathcal{S}}$. Here, $\lambda$ is a control parameter, $\xi^{i}_{q}\in[-1,1]$ is an arbitrary real number, and $X^{i}_{q}$ is the $i^{th}$ single-qubit Pauli operator for the qubit $q$ ($1=X$, $2=Y$, and $3=Z$). We assume all $X_{q}^{i}$ anti-commute with at least one member of $\mathcal{S}$. We note that one could replace $X_{q}^{i}$ with any anti-commuting set of Pauli operators with local support (we are implicitly defining local as any Pauli operator with support less than that of any stabilizer). We pick all single-qubit operators for simplicity, although our results generalise straightforwardly. We also note that a large class of stabilizer code Hamiltonians are CSS codes, for which the Hamiltonian schematically takes the form

[TABLE]

where $O_{s}^{X}$ ($O_{s}^{Z}$) is a product of only $X$-type ($Z$-type) Pauli operators. While all models considered in this paper are CSS codes, our results do not rely on this assumption. In what follows, we first show that in the absence of any perturbations, the entanglement spectrum of a stabilizer code Hamiltonian $H_{\mathcal{S}}$ is flat, and then proceed to study non-trivial universal features which appear in the ES by considering $V$ as a perturbation to $H_{S}$.

II.2 Eigenstates and Entanglement Spectrum of Unperturbed Stabilizer Codes

Let us first consider a stabilizer Hamiltonian in the absence of any perturbations. As mentioned in the previous section, we can extend $G\to\tilde{G}$ in Eq. (4), for $\tilde{G}$ generated by $\{O_{i}\}_{i\leq d_{G}}$ along with a mutually commuting set of logical operators. Then, $\mathbf{k}\in\{0,1\}^{N}$ and

[TABLE]

i.e. $\mathcal{P}^{\mathbf{k}}$ is the pure state density matrix. For the energy of the state $\ket{\mathbf{k}}$, when $|\mathcal{S}|\neq N$, then the terms $O_{i}$ from our generating set for $G$ in the Hamiltonian each contribute $-(-1)^{k_{i}}$, leading to a total energy

[TABLE]

where $\|\|$ is the Hamming weight. The contributions from the remaining terms ($\mathcal{S}-\{O_{i}\}_{i\leq d_{G}}$) depend on the specific generating set chosen but in general, these contributions will cancel in the energy denominators used in the perturbative expansion (see Sec. II.3).

We now calculate the entanglement spectrum of any eigenstate of the stabilizer Hamiltonian for a bipartition $(A,B)$, where $A$ is some subset of the qubits forming $\mathcal{H}$ and $B=A^{c}$. Following Haah (2013); Ma et al. (2018b); Schmitz et al. (2018), one can show that for any $(A,B)$ such that $A$ is smaller than the code distance333The code distance is the minimum size of the support over all members of $L$., the reduced density matrix for any state $\mathbf{k}$ is

[TABLE]

where $s_{A}=|A|-d_{G_{A}}$ and $d_{G_{A}}$ is the dimension of the subgroup $G_{A}\subseteq G$ which only has support in $A$. $\mathcal{P}_{A}^{\mathbf{k}_{A}}$ is the analogous projection operator to Eq. (4) for $G_{A}$ and $\mathbf{k}_{A}$ is the part of $\mathbf{k}$ which corresponds to $G_{A}$.

We have thus shown that the reduced density matrix is proportional to a projection operator, thereby proving that the ES of any unperturbed stabilizer code is flat i.e., eigenvalues of the reduced density matrix are all equal; in the process, we have also obtained the von Neumann entropy for $A$:

[TABLE]

The same calculation can be carried out from the perspective of $B$ and we would find a nearly identical form of the reduced density matrix $\rho_{B}$ for the analogous $G_{B}$. Correspondingly, one can show that the entanglement entropy for $B$ is

[TABLE]

These results also serve to highlight the relative simplicity and power of the stabilizer formalism, as we have made no assumptions besides a stabilizer description for the system of interest in the derivation. In order to extract interesting physics, one needs to further examine the behaviour of Eqs. (11) and (12) in the presence of perturbations/deformations: for topologically ordered states, one is typically interested in extracting the sub-leading corrections to Eq. (12) which are invariant under arbitrary deformations of the partition, while for SPT states, one expects a non-trivial degeneracy in the ES which is robust against arbitrary local perturbations respecting the symmetry protecting the bulk state. We relegate the discussion of such non-trivial features to Secs. III and IV, where we analyse the ES of various models in the presence of perturbations.

A related concept to the entanglement entropy for stabilizer codes is that of the recoverable information, defined for any stabilizer code and bipartition $(A,B)$ as

[TABLE]

where $d_{\partial}$ is the number of members of $\{O_{i}\}_{i\leq d_{G}}$ which are “cut,” i.e., have support in both $A$ and $B$, and the minimization is over all possible choices of a generating set, assuming open boundary conditions. In Ref. Schmitz et al. (2018), it was shown that $\mu>0$ and, generally, that $d_{\partial}-s_{A}-s_{B}$ is equal to the dimension of the non-local surface stabilizer (NLSS) group $G_{\text{NLSS}}$ defined as

[TABLE]

where $G_{\partial}$ is the group generated by the cut members of $\{O_{i}\}_{i\leq d_{G}}$. Thus, the recoverable information is the minimum dimension over all NLSS groups. Every member of $g_{\partial}\in G_{\text{NLSS}}$, referred to as an NLSS, has the general form

[TABLE]

where $g_{A}\in G_{A}$ and $g_{B}\in G_{B}$. We can interpret each of these as an emergent Gauss’ law constraint by noting that if $g_{A}=\prod_{i\in F_{A}}O^{A}_{i}$, for stabilizers $O_{i}^{A}\in G_{A}$ and $F_{A}$ a subset of stabilizer indices in $A$, then by restricting the NLSS to its support in $A$ (where we note $(g_{B})_{A}=I_{A}$) we have

[TABLE]

Thus, a measurement purely on the boundary is equal to the number of “charges” in a subset of the bulk (mod 2). In general, either $g_{A}$ or $g_{B}$ has the given form, but it may happen that $g_{A}=I_{A}$ and $g_{B}\neq I_{B}$, or vice versa, i.e.,

[TABLE]

or some product of cut stabilizers is only supported in $B$. We refer to such an NLSS as a superficial NLSS since the corresponding Gauss’ law is only along the boundary. Every constraint which contains a cut stabilizer always implies an NLSS. However, not all NLSS are formed this way. The minimization in the definition of recoverable information removes all NLSS coming from trivial constraints with the aim of capturing only those arising from topological constraints. We return to the importance of NLSS and how they protect the emergent Gauss’ laws in the ES in Sec. V.

Before proceeding to the analysis of the ES of stabilizer codes in the presence of perturbations, we need to establish some further properties of the eigenstates $\ket{\mathbf{k}}$, for which we can infer a Schmidt decomposition for a given cut $(A,B)$ from Eq. (11). Any stabilizer in $A$, $O^{A}_{i}\in G_{A}$, maintains the same eigenstate with respect to the reduced density matrix $\rho_{A}$, thus allowing us to partially index all Schmidt vectors by their eigenvalues for the stabilizers in $A$. We then require a basis for the subspace defined by the projectors $\mathcal{P}_{A}^{\mathbf{k}_{A}}$, and we use cut stabilizers for this purpose. Likewise, we do the same for $G_{B}$

Since all of $G_{A}$ is supported in $A$, the support of any cut stabilizer in $A$ necessarily commutes with all of $G_{A}$. Further, as the recoverable information is always positive, there exist more cut stabilizers than $\log_{2}$ of the rank of $\mathcal{P}^{\mathbf{k}_{A}}_{A}$. Thus, it is always possible to choose a mutually commuting set of cut stabilizers (with no unique choice for this set) whose simultaneous eigenstates span this space and (nearly) suffice to completely label the Schmidt vectors. Let $G^{R}_{\partial}\subseteq G_{\partial}$ be the reduced boundary group generated by the chosen basis for the cut stabilizers. Members of $G^{R}_{\partial}$ can be indexed by $\mathbf{k}_{\partial}\in\{0,1\}^{s}$, where $s=s_{A}=s_{B}$ is the entanglement entropy. Then the Schmidt decomposition can be written as

[TABLE]

with

[TABLE]

where $O_{i}$ is a generator of $G^{R}_{\partial}$, and likewise for $B$. Note that the Schmidt decomposition Eq. (19) is not unique as a result of the flat ES; the choice made here is informed by the fact that $\ket{\mathbf{k}}$ must have the same stabilizer eigenvalues in its Schmidt form. This is why the two boundary indices in each term must have a binary sum of $\mathbf{k}_{\partial}$. Likewise, $\mathbf{M}$ is an $s\times d_{G}$ binary matrix which maintains the eigenvalues for the stabilizers in $G_{\partial}-G^{R}_{\partial}$. Such operators must anti-commute with both pieces of some cut stabilizers in $G_{\partial}$ (since it commutes with the complete stabilizers). This flips the $\mathbf{l}_{\partial}$ eigenstates by some $\mathbf{p}^{i}\in\{0,1\}^{s}$. Upon relabelling the sum, we find that $\mathbf{M}$ must satisfy $\mathbf{p}^{i}\mathbf{M}\mathbf{k}=\mathbf{k}_{i}$.

II.3 Perturbed Density Matrix and Entanglement Hamiltonian

Thus far, we have shown that, in the absence of any perturbations, the reduced density matrix for a stabilizer Hamiltonian (2) is proportional to a projection operator; consequently, the unperturbed entanglement spectrum is flat. We now consider the Hamiltonian (6), which includes perturbations of the form specified in Eq. (7), to identify universal features of the ES which persist in the presence of such perturbations. Before delving into the full derivation, we outline the general idea underlying our method for finding an approximate density matrix for the perturbed ground state i.e., the ground state of Eq. (6). We make the following general assumptions:

• •

The perturbation is weak enough so as to not close the energy gap i.e., we remain in the same phase described by the exactly solvable stabilizer Hamiltonian.

• •

We work in the thermodynamic limit of the system.

• •

We require that the sub-region $A$, although much smaller than the entirety of the system, is much larger than the size of the stabilizers (region over which a stabilizer has non-trivial support) i.e., the linear size $R$ of region $A$ is assumed to obey $\xi\ll R\ll L$, where $\xi$ is the correlation length of the perturbed Hamiltonian and $L$ is the linear size of the system.

We make use of unitary perturbation theory (see Appendix A), which is a variation on the oft-utilised Schrieffer-Wolff perturbation theory Schrieffer and Wolff (1966); Bravyi et al. (2011). The central focus of unitary perturbation theory (UPT) is to approximate a unitary operator $U$ which maps unperturbed eigenvectors $\ket{n}$ to eigenvectors $\ket{n^{\prime}}$ of the perturbed system. As with any controlled perturbative scheme, this can be done up to some fixed order in the perturbative parameter $\lambda$. Where UPT differs from conventional perturbation theory is that while the latter truncates the expansion of the system state $\ket{\psi}$ to a given order in $\lambda$, the former instead truncates the expansion of the anti-hermitian generator of the unitary $U$ to that order in $\lambda$. This ensures that the approximate transformation $U$ maintains unitarity at any finite order in its perturbative expansion. To first order in the control parameter $\lambda$, the unitary operator is given by $U=\exp(\lambda\mathcal{L})$, with

[TABLE]

and where $V_{\mathbf{k}\mathbf{l}}=\braket{\mathbf{k}}{V}{\mathbf{l}}$ and $[.]$ is the Iverson bracket which equals 1 if the proposition inside is true and 0 otherwise (see Appendix A). Making a further approximation valid in the large $A$ limit, we consider the resulting action of the unitary $U$ on the Schmidt form Eq. (19), allowing us to perform the partial trace required to form the perturbed density matrix $\tilde{\rho}_{A}$ for the ground state of the perturbed Hamiltonian. The perturbed entanglement Hamiltonian (EH), to first order, is then defined as

[TABLE]

where $\simeq$ implies unitary equivalence. Hence, $\lambda$ may be thought of as the inverse temperature for this state.

We now describe the derivation of $H_{\text{ent}}$ in detail, starting with the matrix elements of the perturbation, Eq. (7). Recall that $X_{q}^{i}$ is the $i^{th}$ single-qubit Pauli operator for qubit $q$, with $1=X,2=Y,$ and $3=Z$. For every $X_{q}^{i}$, we can assign a binary string $\mathbf{p}_{q}^{i}\in\{0,1\}^{N}$ representing all stabilizer basis elements which anti-commute with $X_{q}^{i}$. Thus,

[TABLE]

where $\oplus$ is the binary sum or bitwise-XOR of the two strings. From this we get the modulus-square of each perturbation term

[TABLE]

However, this does not resolve the phase of the matrix element: $X_{q}^{i}\ket{\mathbf{k}}\propto\ket{\mathbf{k}\oplus\mathbf{p}_{q}^{i}}$, up to an overall phase. Nonetheless, one can complete a canonical basis for all Pauli operators as described in Appendix B using the stabilizer (and logical) basis operators and some canonical duals. This allows us to write all $X_{q}^{i}$ in terms of this basis, which implies the existence of a $\tilde{\mathbf{p}}_{q}^{i}\in\{0,1\}^{N}$ such that

[TABLE]

Thus, the matrix element $V_{\mathbf{k}\mathbf{l}}$ is given by

[TABLE]

We can also rewrite the energy denominator as

[TABLE]

where the overbar denotes the complement, or negation, of the string.

Putting the terms together, we see that the generator of the perturbation unitary is

[TABLE]

where we have defined the $C$-coefficient

[TABLE]

Despite the apparent complexity of the $C$-coefficients, a crucial property which can be easily established is that for qubits $q,r$ which are not both contained within the support of some stabilizer,

[TABLE]

which reflects the fact that if two qubits are far enough separated, they do not interact. This becomes clear upon examination of the energy denominator, which is sensitive only to the string in the vicinity of $q$. In other words, the binary dot product projects $\mathbf{k}$ and $\mathbf{k}\oplus\mathbf{p}_{r}^{j}$ onto $\mathbf{p}_{q}^{i}$ and, if the two qubits $q,r$ are not both contained within the support of some stabilizer (are far enough away), the projection is not affected by the change $\mathbf{k}\to\mathbf{k}\oplus\mathbf{p}_{r}^{j}$. As for the overall phase, the invariance of the $C$-coefficients under this change can be seen from

[TABLE]

where we have used the (anti) commutation relations for the strings under the assumption $q\neq r$. Another important property of the $C$-coefficients can be established from the requirement that $\mathcal{L}$ is skew-Hermitian:

[TABLE]

Expressing the skew-Hermitian generator of the perturbation unitary in terms of the $C$-coefficients and using their aforementioned properties allows us to expand the perturbation unitary $U$ as

[TABLE]

where, for ease of notation, we have suppressed the Pauli-type index, and defined $\mathbf{p}_{\vec{q}}\equiv\oplus_{i}\mathbf{p}_{q_{i}}$ and

[TABLE]

Eq. (33) can be verified through induction.

In order to perform the partial trace, we make a further approximation that for any string $\vec{q}$, we can reasonably make the substitution,

[TABLE]

where $\vec{q}_{A(B)}$ is the part of the string $\vec{q}$ contained in region $A(B)$ while maintaining their relative order. While the binary sum in the bra, $\mathbf{p}_{\vec{q}}=\oplus_{i}\mathbf{p}_{q_{i}}=\mathbf{p}_{\vec{q}_{A}}\oplus\mathbf{p}_{\vec{q}_{B}}$ is exact, the corresponding coefficient in the expansion of the perturbation unitary $U$ is not, making this a non-trivial approximation whose effect must be closely evaluated. From the definition of the $C$-coefficients Eq. (29), the fact that $\xi_{q}\in[-1,1]$, and that the energy gap is always bounded by $1$, it is straightforward to see that

[TABLE]

Hence, to evaluate the effect of the approximation Eq. (35), we must consider how many terms in the sum are in error. It should be clear from the definition of the $C$-coefficients that only the energy denominators need to be considered; since the energy denominators “project” on the $\mathbf{k}$ string in the vicinity of $q$, if $q_{1}$ and $q_{2}$ are not shared within the support of any stabilizer, Eq. (30) implies that $C^{q_{2}}_{\mathbf{k}\oplus\mathbf{p}_{q_{1}}}=C^{q_{2}}_{\mathbf{k}}$. At first order, our approximation is clearly exact. At second order, suppose that we are considering a $d$-dimensional system such that the size of the system $\sim L^{d}$ for $L$ large and that the size of region $A\sim R^{d}$ with $R<L$. Of the $\sim L^{2d}$ total terms, only $\sim R^{d-1}$ are in error as only this number of terms have both qubits in the boundary and in the same cut stabilizer. The number of terms in error thus comprise a vanishing, negligible fraction in the thermodynamic limit. On the other hand, one might worry that only terms which involve the boundary will matter. However, even considering only those terms for which both qubits are in the boundary, the ratio of error terms to non-error terms goes as $R^{-(d-1)}$, and so the error is suppressed as $R$ becomes large (assuming $d\geq 2$). The suppression of higher-order errors follows by a similar logic.

Effectively, the replacement (35) amounts to the assumption that all terms for which $\vec{q}$ is composed of $\vec{q}_{A}$ and $\vec{q}_{B}$ are equal. We must thus account for the multiplicity of these terms. For $\vec{q}_{A}$ of length $\beta$ and $\vec{q}_{B}$ of length $\alpha-\beta$, the multiplicity is $\binom{\alpha}{\beta}$.

We now apply the perturbation unitary $U$ to a ground state $\ket{0}$ of the unperturbed Hamiltonian in order to find the approximate ground state for the perturbed Hamiltonian (6). Using the Schmidt form of the unperturbed eigenvectors Eq. (19), and within the approximation (35), we find that

[TABLE]

where by $\left(\mathbf{p}_{\vec{q}_{A}}\right)_{A\cup\partial A}$ we mean the the part of $p_{\vec{q}_{A}}$ corresponding to $G_{A}$ and $G_{\partial}^{R}$, and likewise for $B$. $\mathbf{k}_{\partial}$ is short for $(\mathbf{0}_{A},\mathbf{k}_{\partial})$, which labels members of $G_{A}\oplus G^{R}_{\partial}$. Note in particular that we have split-up the binary values of $\mathbf{p}_{\vec{q}_{A}}\oplus\mathbf{p}_{\vec{q_{B}}}$ such that if $q\in A$, then it only contributes to $G_{A}\oplus G^{R}_{\partial}$ in the Hilbert space of $A$. This is always possible since an edge in $A$ can only anti-commute with members of $G_{A}\oplus G^{R}_{\partial}$. Likewise is true for $B$. This splitting ensures that the binary sum of the values in the $A$ and $B$ kets equals the value in the original ket, as per our discussion in Section II.2.

While the final expression obtained in Eq. (II.3) is ostensibly another Schmidt form for the approximate perturbed ground state, this assumes that $\{\ket{\psi_{\mathbf{k}_{\partial}}}_{A}\}$ and $\{\ket{\phi_{\mathbf{k}_{\partial}}}_{B}\}$ are orthonormal. Focusing on $A$, we consider the overlap of any two of these states,

[TABLE]

If $\mathbf{k}_{\partial}=\mathbf{l}_{\partial}$, then $\left(\mathbf{p}_{\vec{q}_{1}}\right)_{A\cup\partial A}=\left(\mathbf{p}_{\vec{q}_{2}}\right)_{A\cup\partial A}$ within the sum, causing the phase factors to cancel and all dependence on $\mathbf{k}_{\partial}$ drops out. This implies that $\braket{\psi_{\mathbf{k}_{\partial}}}{\psi_{\mathbf{k}_{\partial}}}$ is a constant which, without loss of generality, we take to be $1$. However, when $k_{\partial}\neq l_{\partial}$, the fact that the $C$-coefficients are insensitive to $\mathbf{k}_{\partial},\mathbf{l}_{\partial}$ implies that the inner product does not sum to zero. Nonetheless, we can still make use of this form to perform the partial trace, which is given by

[TABLE]

where $\mathcal{Z}$ is the normalisation factor which makes $\tilde{\rho}_{A}$ have unital trace after the approximation is made. To further simplify the expression, we define

[TABLE]

and analogously for $\tilde{U}_{B}$. Using these definitions, we can write

[TABLE]

and

[TABLE]

where

[TABLE]

with $\mathcal{P}_{A}^{\mathbf{0}}$ ($\mathcal{P}_{B}^{\mathbf{0}}$) a projector for $G_{A}(G_{B})$ (see Sec. II.2).

Note that unitary equivalence in Eq. (II.3) is a result of the fact that, for two operators $O_{1}$ and $O_{2}$, $O_{1}O_{2}$ has the same non-zero eigenvalues as $O_{2}O_{1}$. Since we are only interested in the density matrix up to first order in the control parameter $\lambda$ and since $\sigma_{B}\in\mathcal{O}(\lambda)$, we can expand the square root as $\sqrt{\mathcal{P}_{A}^{\mathbf{0}}+\sigma_{B}}=\mathcal{P}_{A}^{\mathbf{0}}+\frac{1}{2}\sigma_{B}+\mathcal{O}(\lambda^{2})$. Again, to lowest order we find that

[TABLE]

where we have made use of the fact that $\tilde{U}_{A}$ is the identity at lowest order and that $\mathcal{P}_{A}^{\mathbf{0}}\sigma_{B}=\sigma_{B}\mathcal{P}_{A}^{\mathbf{0}}=\sigma_{B}$.

We are now well positioned to define the two parts of the EH:

[TABLE]

The $A,B$ symmetry manifest in these two terms should be unsurprising given that the entanglement spectrum is the same for both $A$ and $B$, with contributions naturally arising from both sides of the boundary. This also suggests that the non-triviality of the EH is a consequence of the non-unitarity of the operators $\tilde{U}_{A}$ and $\tilde{U}_{B}$. Evaluating the expression (45a), we find in general that

[TABLE]

from which we immediately find that

[TABLE]

Equations (II.3) and (II.3) constitute the central results of this paper, upon which we will now elaborate.

Let us consider $H_{\text{ent}}^{A}$ (II.3), where for $q\in A$ $\tilde{X}_{q}=\mathcal{P}_{A}^{\mathbf{0}}X_{q}\mathcal{P}_{A}^{\mathbf{0}}$, such that the only perturbations which survive are those which commute with all of $G_{A}$—hence, only members of $\partial A$ need be considered. Crucially, we note that there often exist constraints on the terms of the (perturbative) EH: if any product of terms forms a member of $G_{A}$, i.e., there exists a subset $F_{\partial}\subseteq\partial A$ such that $\prod_{q\in F_{\partial}}X_{q}\in G_{A}$, the projection must satisfy

[TABLE]

This constraint is a consequence of the NLSS, as discussed in Sec. II.2 where $\prod_{q\in F_{\partial}}\tilde{X}_{q}=\left(g_{\partial}\right)_{A}=\prod_{i\in F_{A}}\mathcal{O}_{i}$ for some set of indices $F_{A}$, as described in Eq. 17. Such a constraint indicates the presence of non-trivial entanglement features in the model under consideration. Equivalently, the constraint (48) may be thought of as defining a $\mathbb{Z}_{2}$ topological surface charge, with the entanglement Hamiltonian confined to the zero charge sector; we return to a closer examination of these features in Sec. V.

Let us now turn to $H_{\text{ent}}^{B}$ (II.3), for which we must consider the form of the operators $\tilde{Z}_{q}$ for $q\in B$. If $X_{q}$ commutes with all of $G_{B}$, then its action must flip only cut stabilizers, implying that $\tilde{Z}_{q}$ is an operator which flips the same cut stabilizers, but with respect to their support in $A$. Such an operator can always be found by considering a canonical basis for Pauli operators in $A$ by using a set of $|A|$ complete and cut stabilizers in $A$, and then finding their canonical duals. Thus, $\tilde{Z}_{q}$ is necessarily the Pauli operator formed as the product of all operators dual to the cut stabilizers with which $\tilde{Z}_{q}$ anti-commutes. As constructing $\tilde{Z}_{q}$ in this manner tends to be a fairly tedious procedure in practice, we now discuss a simpler, more efficient method.

Since $\tilde{Z}_{q}$ is a Pauli operator, let us consider $Z_{q}X_{q}$ (where the lack of a tilde signifies that we are considering the corresponding Pauli operator in the full space and without the projection). As $Z_{q}$ and $X_{q}$ anti-commute with the same members of $\mathcal{S}$, their product necessarily commutes with all of $G$. This implies $Z_{q}X_{q}\in\tilde{G}$, i.e., it is either a product of stabilizers or a logical operator. If $A$ is significantly smaller than the code distance and $X_{q}$ is local, we can generally conclude that $Z_{q}X_{q}$ can not be a logical operator. Thus $Z_{q}X_{q}$ must be a member of $G_{\partial}$ such that its only support in $B$ is $X_{q}$ and only support in $A$ is $Z_{q}$ i.e. it represents a cut stabilizer group element. Further, just as was the case with $\tilde{X}_{q}$, all $\tilde{Z}_{q}$ are subject to the same constraint that some product of these operators must act as $\mathcal{P}_{A}^{\mathbf{0}}$ if they form an NLSS. By definition, if there exists a set $F_{\partial}\subseteq\partial B$ such that $\prod_{q\in F_{\partial}}X_{q}$ forms an NLSS in $G_{B}$, $\prod_{q\in F_{\partial}}Z_{q}$ forms the corresponding NLSS in $G_{A}$ as every NLSS has the form $g_{\partial}=g_{A}g_{B}$ and $g_{A}$ is unique to $g_{B}$ (see Ref. Schmitz et al. (2018)). Note that the only guarantee that $Z_{q}$ is local is if $X_{q}Z_{q}\in\mathcal{S}$, i.e., the product forms a single stabilizer.

We note that even though we included only single-qubit Pauli operators in the perturbation $V$ to the stabilizer Hamiltonian, if all such terms of a given Pauli type do not survive the projection, we are then forced to consider higher order perturbations. Rather than computing second-order corrections coming from UPT, we can instead add local perturbations consisting of two-qubit Pauli operators. All of the results from above follow in a similar fashion, with the caveat that we now have to carefully consider those perturbation terms which are themselves cut. However, since every single-qubit Pauli does not survive, such cut perturbation terms will also not survive and can hence be safely ignored. We refer to such a process as second-order even though it technically arises at first-order in UPT, albeit from two-qubit perturbations444The use of this term is justified as, even for single-qubit Pauli operator perturbations, we expect that the EH contains higher-order terms of the same form but with different coefficients. In other words, second order UPT with single qubit Pauli operator perturbations yields the same results, up to unimportant coefficients, as first order UPT with two qubit Pauli operator perturbations. Our choice to proceed with the latter is made simply to avoid a lengthy digression into the derivation of second order UPT coefficients.. Analogously, if all second-order contributions fail to survive, we continue to third-order contributions and so forth, until all lowest order contributions are identified.

II.4 Entanglement Hamiltonian for different topological sectors

Until this point, we have effectively assumed the ground state under consideration is the $+1$ eigenstate under all logical operators. For $A$ smaller than the code distance, the cleaning lemma Bravyi and Terhal (2009); Haah (2016) ensures that no logical operator is cut and thus its eigenstate does not affect the resulting entanglement Hamiltonian. However, if $A$ is a non-cleanable subset of the qubits, then some logical operators are necessarily cut. What this means is simply that some subset of the terms may not form NLSS’s but rather form logical operators. If the original ground state is in the $-1$ eigenstate of some logical operator, then this has the effect of projecting the entanglement Hamiltonian onto the one-charge sector for the topological surface charge defined by that logical operator. Since both NLSS’s and logical operators can be connected to topological constraints amongst the stabilizers (see Refs. Schmitz et al. (2018); Schmitz (2018)), we can define these topological surface charges in a bipartition-independent way (i.e., a topological surface charge is defined by its relation to a topological constraint). This also has the rather interesting implication that a topological surface charge can always measure a topological charge in the bulk, which constitutes a bulk-boundary correspondence for topologically ordered phases. We return to a discussion of topological surface charges and their relation to topological constraints in Sec. V.

III Entanglement Spectra for Long-Range Entangled States

Having established the general formalism for deriving the entanglement spectra for stabilizer codes in the presence of generic perturbations, we now analyse the resultant entanglement Hamiltonians given by Eqs. (II.3) and (II.3) for specific models. In this section, we discuss stabilizer codes describing LRE quantum phases of matter, which include the $d=2$ toric code/Wen-plaquette model and $d\geq 3$ toric code as examples of conventional topological order. We then consider the X-cube model Vijay et al. (2016) and Haah’s cubic code Haah (2011) as examples of type-I and type-II fracton orders respectively555In type-I models, fractons are created at the corners of membrane-like operators, while they are created at the corners of fractal-like operators in type-II models. Type-II phases are also distinct from their type-I counterparts by the absence of any topologically non-trivial mobile quasi-particles, while type-I models generally host sub-dimensional particles in addition to fractons.. We discuss the existence of an edge-ES correspondence for these systems towards the end of this section.

Although our general derivation of the EH in the previous section did not rely upon the simplifying assumption of considering only CSS codes, all models studied in this section are CSS stabilizer Hamiltonians of the form Eq. (8). We also note that while we describe the NLSS below, a more thorough discussion for most cases considered here can be found in Ref. Schmitz et al. (2018). As a matter of notation, if a stabilizer is cut, we refer to it as an $(a,b)$ cut stabilizer, where $a(b)$ is the size of its support in $A(B)$ for a bipartition $(A,B)$ of the qubits forming $\mathcal{H}$, the Hilbert space of the system in question. Throughout this paper, we will only consider bipartitions associated with the degrees of freedom living in spatially distinct regions $(A,B)$ (with $B=A^{c}$), such that the boundary between regions $A$ and $B$ defines the entanglement cut.

III.1 $\mathbb{Z}_{2}$ Topological Order

III.1.1 $d=2$ Toric Code/Wen-plaquette Model

We start with the $d=2$ toric code Kitaev (2003) and the Wen-plaquette model Wen (2003), examples of phases hosting $\mathbb{Z}_{2}$ topological order which are well-known to be unitarily equivalent in the bulk via a local unitary transformation and a $\pi/4$ rotation of the lattice (see e.g. the discussion in Ref. Wen (2004)). Hence, we only discuss the toric code explicitly here, with results for the Wen-plaquette model following immediately. The toric code in $d=2$ spatial dimensions is defined on the square lattice with one qubit associated with each edge and is described by the Hamiltonian

[TABLE]

The first term in the Hamiltonian is associated to every vertex $v$ of the lattice such that $A_{v}=\prod_{e@v}X_{e}$, i.e., the $X$-type Pauli for each edge attached to $v$. The second term is associated with every plaquette $p$ such that $B_{p}=\prod_{e\in p}Z_{e}$, i.e., every $Z$-type Pauli forming $p$.

For generic entanglement cuts, any boundary666To avoid verbosity, and since it should be clear from context, we will refer to a “boundary defined by the entanglement cut” as simply a “boundary,” except when contrasting it with a physical boundary of the system. parallel with the coordinate directions contributes to the EH at first order, since both vertex and plaquette stabilizers are either $(3,1)$ or $(1,3)$ cut stabilizers. Analogously, any corners of the entanglement cut contribute at second order as both stabilizer types are $(2,2)$ cut stabilizers; this is also the case for straight boundaries which cut diagonally with respect to the coordinate directions (which is the natural cut for the Wen-plaquette model).

We now consider these two types of local boundaries in detail, both of which are depicted in Fig. 1. The first is the case of a locally flat edge which contributes terms to the EH of the form,

[TABLE]

where $e_{A}\in A$ ranges over all edges parallel to the boundary $\partial$ of the entanglement cut and $e_{B}\in B$ ranges over all perpendicular edges. $v(e_{B})$ corresponds to the vertex in $A$ which is attached to the edge $e_{B}$ and $\tilde{T}_{v}=\left(A_{v}\right)_{A}$ is the T-shaped cut vertex stabilizer at $v\in\partial$ (see Fig. 1). The second case is along a “stair-stepping” edge; since there are no contributions at first order, to second order we find that the EH contains terms of the form

[TABLE]

where $v$ ranges over all vertices on the boundary, and $p(v)$ is the plaquette such its intersection with the boundary is $v$. Here, $\tilde{E}_{v}^{X}=\left(A_{v}\right)_{A}$ and $\tilde{E}_{p}^{Z}=\left(B_{p}\right)_{A}$ are the “elbow”–shaped operators formed by both the vertex and plaquette operators (see Fig. 1). Considering these two possibilities is sufficient to obtain all contributions to the EH coming from any bounded entanglement cut.

Considering the commutation relation between the various terms, we see there exists a mapping from the resulting EH to the $d=1$ transverse-field Ising model (TFIM). In particular, any $X$-type operator (be that an edge or a concave elbow operator) anti-commutes with the two neighboring $Z$-type operators (be they T-shaped or convex elbow operators) and vice versa. These are precisely the commutation relations between the terms of the TFIM. This implies that we can map the spin-$1/2$ degrees of freedom along the boundary of any region $A$ to effective spin-$1/2$ degrees of freedom along a $d=1$ cycle of the same size as the boundary, whereby the EH is mapped onto a TFIM:

[TABLE]

where $\tilde{\mathcal{Z}_{i}},\tilde{\mathcal{X}_{i}}$ are the Pauli operators acting on the effective spin degrees of freedom, $\xi_{i},\xi_{ij}$ are determined by the exact form of the mapping, and the tilde represents that we still have to consider the projection due to the presence of NLSS.

For each stabilizer type, there is an NLSS guaranteed by a topological constraint, i.e. the product of all stabilizers of that type is the identity. Each NLSS takes the appearance of a closed string wrapping the boundary and is given by the product of all stabilizers of that type inside of $A$. These NLSS constrain the terms of the EH such that

[TABLE]

with the projectors $\mathcal{P}^{0}_{A}$ defined in Sec. II.3. This implies that after the mapping onto the TFIM, the effective spins are constrained by

[TABLE]

The above constraints imply that the NLSS “enforce” a $\mathbb{Z}_{2}$ invariance, or charge conservation, on the EH. Note that the well-known Kramers-Wannier, or self duality of the $d=1$ TFIM, which maps the transverse field terms $\mathcal{X}$ onto the Ising terms $\mathcal{ZZ}$, is exact only up to the global constraint, i.e. Eq. (54a) is automatically enforced whereas Eq. (54b) is not. However for the EH (52), the duality is exact since both terms of the EH satisfy the constraints as enforced by the NLSS. There is thus an ambiguity in the mapping of the terms which can be understood as a consequence of the $e-m$ duality in the bulk.

From the preceding discussion, we see that the EH for the $d=2$ toric code shows universal behaviour, since, for generic cuts and for arbitrary weak local perturbations, it can be mapped onto a $\mathbb{Z}_{2}$ invariant TFIM of the form Eq. (52), which acts on a chain of effective spin-$1/2$ degrees of freedom. Importantly, it is the NLSS which enforce the global $\mathbb{Z}_{2}$ invariance of the EH, which can thus be understood as a consequence of the bulk $\mathbb{Z}_{2}$ topological order. Similarly, we see from the NLSS that the electromagnetic duality of the bulk state manifests itself as Kramers-Wannier duality in the EH, which suggests that there exists a finite region in parameter space for which the EH can be mapped onto the critical Ising CFT. In particular, for entanglement cuts corresponding to a flat boundary, imposing translation symmetry along the boundary forces all $\xi$ coefficients in the perturbation $V$ to be equal, which maps the EH onto the critical 1d Ising CFT. The condition that all $\xi$ coefficients be equal for the stair-stepping edge can be understood as being enforced by translation symmetry along a flat entanglement boundary of the Wen-plaquette model, in which case the EH again maps onto the 1d Ising model at its critical point. Since this is not expected to hold for generic perturbations, we emphasise that the key universal property of the EH is its mapping onto a $\mathbb{Z}_{2}$ invariant Hamiltonian acting on effective spin-$1/2$ degrees of freedom, with the $\mathbb{Z}_{2}$ invariance enfored by the bulk topological order through the NLSS. We note that our results for the EH of the toric code/Wen-plaquette model are in agreement with those found in Ref. Ho et al. (2015).

Finally, for the case where $A$ is not bounded, but rather spans one direction of the torus, all results remain the same except for the fact that the NLSS get promoted to logical operators. As a consequence, the right-hand sides of Eqs. (53) are replaced by $(-1)^{\ell_{1}}(-1)^{\ell_{2}}\mathcal{P}_{A}^{(\ell_{1},\ell_{2})}$, where $\ell_{1},\ell_{2}$ are the indices for the logical operators. When $A$ was bounded, the projection was necessarily onto the “zero charge” sector of the TFIM, but due to the possible negative sign, the projection is now onto the charge sector which corresponds to the topological sector in the bulk. As a consequence, if we are in the 1-charge sector for both topological indices, the $\mathcal{Z}\mathcal{Z}$ terms must multiply to $-1$, something which is otherwise impossible in the 1d TFIM. To understand this better, consider a flat cut extending along one direction. Even though the product of all $T$ operators has no support on the qubits along the boundary (the hypothetical TFIM degrees of freedom), this operator extends into the bulk and can hence have an eigenstate which differs from one.

III.1.2 $d\geq 3$ Toric Code

The most natural, albeit by no means the only, generalisation of the toric code for $d\geq 3$ is defined on the $d$-dimensional Euclidean lattice, with one qubit per edge. The Hamiltonian $H_{TC_{d}}$ is identical to that of the $d=2$ toric code (49), with the stabilizers of the exact same form, where the vertex stabilizer $A_{v}$ is formed by the $2d$ X-type Pauli operators attached to a vertex and the plaquette stabilizer $B_{p}$ is still the product of the four Z-type operators forming a plaquette.

We focus here on the $d=3$ toric code and begin by looking at contributions to the EH stemming from the vicinity of a flat boundary with respect to the coordinate directions. These parts of the boundary have contributions coming from the $(1,3)$ cut plaquettes which intersect the boundary along an edge as well as from the $(5,1)$ cut vertex stabilizers, as depicted in Fig. 2. These contributions take the form

[TABLE]

Along any “hinges” of the entanglement cut, we have additional $(2,2)$ cut plaquettes for concave hinges and $(4,2)$ cut vertex contributions at convex hinges, as shown in Fig. 2. Likewise, any convex corner contributes one $(3,3)$ cut vertex stabilizer whereas concave corners contribute no additional terms. The presence of such hinges and corners in the entanglement cut add terms to the EH of the form:

[TABLE]

The terms considered above are sufficient for generating all contributions to the $d=3$ toric code EH coming from any bounded entanglement cut. An example of such contributions are shown in Fig. 2 for a cubic entanglement cut.

To understand the universal features of the EH of the $d=3$ toric code, we now consider the NLSS. In general, there exist an extensive number of NLSS since the local cube constraints ($\prod_{p\in\text{cube}}B_{p}=I$) imply that every product of $\tilde{Z}_{e}$ and $\tilde{E}^{Z}_{p(v)}$ which forms a closed loop along the boundary of the cut must equal the identity (see Fig. 2). Note that this includes closed loops which are deformed around hinges and corners. The same constraint is present among the the $\mathcal{Z}\mathcal{Z}$ terms of a $d=2$ TFIM. Similarly to the $d=2$ toric code, this suggests a mapping from the original qubits along the boundary onto effective spin-$1/2$ degrees of freedom on the convex vertices of the boundary 777As a consequence of fixing some degrees of freedom along the boundary—specifically, those associated with the complete plaquette stabilizers—the total Hilbert space dimension of the effective degrees of freedom on the boundary is reduced. This same dimensional reduction is also present in the X-cube model, which we discuss below.. Correspondingly, at the operator level the mapping is implemented through

[TABLE]

where $e=(i,j)$ for $Z_{e}$ terms or $i$ and $j$ corresponds to the two adjacent corners of $p$ for $E^{Z}_{p}$ terms, and $v=i$. As before, $\tilde{\mathcal{Z}}_{i},\tilde{\mathcal{X}}_{i}$ are the Pauli operators acting on the effective spin degrees of freedom and the tilde represents projection due to the presence of NLSS.

The above mapping preserves all the commutation relations between the different terms of the EH and thus, we find that the EH for the $d=3$ toric code is unitarily equivalent to the $d=2$ TFIM:

[TABLE]

where $\braket{i,j}$ now refers to nearest neighbors on the square lattice and $\xi_{i},\xi_{ij}$ are determined by the precise form of the mapping. Note that unlike the $d=2$ toric code, we were able to make the mapping explicit in this case as there is no ambiguity in how the different terms get mapped. This is due to the lack of an exact $e-m$ duality in the $d=3$ toric code, which translates into the absence of an exact self-duality in the EH, as evinced also from its equivalence to the $d=2$ TFIM. While the $d=2$ TFIM does not harbour a Kramers-Wannier duality, it is dual to the $d=2$ Ising gauge theory via the celebrated Wegner duality Wegner (1971), indicating that the EH for the $d=3$ toric code is dual to a theory with a local $\mathbb{Z}_{2}$ symmetry. Further, although the existence of a phase transition in the $d=2$ TFIM shows that there exists a region in parameter space for which the EH is at the critical point, generically we do not expect the EH to be near criticality.

Aside from the “trivial” NLSS, there are also those connected to the topological constraints in the bulk. For topologically trivial entanglement cuts there is only one such constraint, and thus only one NLSS constraint( $\prod_{v}A_{v}=I$) not already captured by the trivial NLSS, namely:

[TABLE]

which for the effective TFIM implies

[TABLE]

This constraint enforces global $\mathbb{Z}_{2}$ invariance on the EH and can be understood as a consequence of charge conservation for the electric sector in the bulk. Much like the EH for the $d=2$ toric code, we hence find that the universal features of the EH for the $d=3$ toric code are its mapping onto a globally $\mathbb{Z}_{2}$ symmetric Hamiltonian acting on effective spin-$1/2$ degrees of freedom, with the symmetry enforced by the bulk topological order through the topologically non-trivial NLSS. Where the EH for the $d=2$ toric code harbours a self-duality, for the $d=3$ toric code the EH is dual to the $d=2$ Ising gauge theory.

So far, it is not clear whether the EH can additionally encode the flux conservation condition in the bulk. In the absence of perturbations, there exist contributions to the entanglement, which appear in the recoverable information, originating from the flux conservation condition. However, this only occurs for topologically non-trivial cuts such that the bulk can enclose or encircle closed flux lines. Since the recoverable information counts non-trivial NLSS, we expect the EH to encode this information as well for topologically non-trivial cuts. We return to this point in Sec. V.

We now consider the case when $A$ is not bounded, i.e., when the boundary wraps around the system. For concreteness, we assume the boundary is a flat surface which is everywhere perpendicular to the $z$-direction. Other possibilities can be analogously understood. All of our discussion for bounded cuts carries over with the key difference being the promotion of non-trivial NLSS to logical operators. The electric charge sector gets projected onto the topological sector given by the quantum number of the membrane logical operator perpendicular to the $z$ direction. Likewise, we also find two flux sectors defined by

[TABLE]

where $\text{path}_{x(y)}$ is any path which wraps around the $x(y)$ direction, and $f_{x(y)}$ is the quantum number for the logical string operator wrapping around the $x(y)$ direction. Similarly to our discussion of the $d=2$ toric code, we find that these conditions are only possible since the mapping of the EH onto the $d=2$ TFIM is exact only up to global constraints.

III.2 Type-I Fracton Order: X-cube model

As the archetypal model displaying type-I fracton order, we now study the X-cube model introduced in Ref. Vijay et al. (2016). The model is defined on a cubic lattice, with qubits living on each edge of the lattice. It is a CSS Hamiltonian formed by two different stabilizer types,

[TABLE]

where $i=x,y,z$. The first stabilizer type $B_{c}$ is associated with each cube $c$ such that $B_{c}=\prod_{e\in c}Z_{e}$ i.e., every Z-type Pauli forming the cube. The second stabilizer type has three stabilizers associated to each vertex, one for each direction, such that $A_{v}^{i}=\prod_{e@v,\perp i}X_{e}$ i.e., every X-type Pauli attached to $v$ and in the plane perpendicular to the direction $i$, as depicted in Fig. 3. Given the extensive literature on fracton order, we do not discuss the properties of this model in detail here but refer the reader to Ref. Nandkishore and Hermele (2018) for a review.

Unlike topologically ordered systems, those with fracton order have a non-trivial geometric dependence Slagle and Kim (2018); Shirley et al. (2018); Prem et al. (2018); Slagle et al. (2018a); Gromov (2018); for the X-cube model, this is evidenced clearly from its sub-extensive ground state degeneracy on the 3-torus Vijay et al. (2016). Due to the geometric sensitivity, tabulating all possible contributions to the EH from the set of all bounded cuts is much more challenging than for the $d=2,3$ toric codes discussed in the previous section. However, such a tabulation is unnecessary since the relevant features distinguishing the entanglement structure of fracton order from that of topological order can be understood by studying an $R\times R\times R$ cubic cut along the coordinate directions, which we now consider.

Within the interior of each plane of this cut, but away from its boundaries (hinges and corners), only those $(3,1)$ vertex stabilizers whose directional index is along that plane survive at first order i.e., only the $A_{v}^{i}$ and $A_{v}^{j}$ $(3,1)$ stabilizers survive on the two $ij$ planes of the cut. We must go to fourth-order before we find a $(4,8)$ cut cube stabilizer inside each plane. Both are demonstrated in Fig. 4. Thus the surface of the cut, away from the hinges, contributes

[TABLE]

where $\tilde{B}_{p}$ is the same plaquette operator as in the toric code and $\epsilon=\lambda^{3}$, since the second term appears at fourth order in perturbation theory. Note that we only include one $\tilde{T}_{v}$ for each vertex along the boundary as the product of the three vertex terms $A_{v}^{i}$ at any given vertex is a constraint. This implies the presence of a trivial NLSS requiring the equivalence of the two cut stabilizers at $v$.

Next, along the hinges of the cut but away from its corners, we find $(1,11)$ cut cube stabilizers which intersect the hinges. Since the cut vertex stabilizers which contribute are the same as those coming from the surface, these have already been accounted for in Eq. (63). In addition, there are also contributions from $(2,2)$ cut vertex (elbow) operators. However, this is a second-order contribution which, due to the superficial NLSS, is the product of two first-order contributions, and is thus excluded. Finally, there are no cube contributions coming from the corners of the cut but there do exist three $(2,2)$ cut vertex stabilizers whose product is the identity. As all three of these appear at the same (second) order, we include all of them. Therefore, we find additional contributions to the EH coming from the hings and the corners, which take the form:

[TABLE]

Eqs. (63) and (64) include all possible contributions to the EH, where the hinge contributions are highlighted in Fig. 4.

Along the surface of the cut but away from hinges and corners, the surface EH, which is formed solely by the terms in Eq. (63), can be better understood by considering the commutation relations between these terms. Every vertex term $\tilde{T}_{v}$ anti-commutes with each of the four plaquette terms $\tilde{B}_{p}$ sharing the vertex $v$; likewise, every $\tilde{B}_{p}$ term anti-commutes with each of the four vertex terms lying on the corners of $p$. We can then define a mapping which takes the original qubits living on edges of the square lattice on the surface onto effective spin-$1/2$ degrees of freedom living on plaquettes of the original lattice, or equivalently, on vertices of the dual square lattice (see Fig. 5). Under this mapping, the $\tilde{T}_{v}$ and $\tilde{B}_{p}$ terms of the surface EH (63) are mapped as follows:

[TABLE]

where $i$ is the plaquette degree of freedom associated with the $p$, and where

[TABLE]

is the product of four (effective) Pauli-$\mathcal{Z}$ operators acting on the vertices $q_{i}$ forming the dual plaquette $q$ (see Fig. 5). As before, the tilde here signify that we must also account for the constraints imposed by the NLSS.

The surface EH, which on the original square lattice along the surface is given by Eq. (63), can hence be mapped onto the subsystem symmetric transverse field Ising-plaquette (TFIP) model

[TABLE]

acting on effective spin-$1/2$ degrees of freedom associated with the dual square lattice on the surface of the cut. On the dual lattice, the $\tilde{\mathcal{B}}_{q}$ term is given by the product of four Pauli-$\mathcal{Z}$ acting on the four vertices forming the plaquette $q$. Unlike the $d=2$ TFIM which has a global $\mathbb{Z}_{2}$ symmetry, the TFIP model instead has a sub-extensive number of $d=1$ subsystem symmetries since it is invariant under flipping all spins along any row or column of the lattice. Thus, the EH for the X-cube model is clearly distinct from that of the $d=3$ toric code.

So far, we have only considered surface contributions to the X-cube EH, but we must also account for the additional hinge and corner contributions given by Eq. (64). While the EH can be mapped onto TFIP models along the six faces of the cut, these models are coupled through additional effective link degrees of freedom living on the “wire-frame” defined by the twelve hinges of the cut. These effective spin-$1/2$ degrees of freedom are then coupled by all terms in Eq. (64) and also by the boundary terms from Eq. (63). Specifically, at the intersection of any two surfaces $s$ and $s^{\prime}$ of the entanglement cut, we defined a mapping such that:

[TABLE]

where $\tilde{\mathcal{B}}_{q}$ was defined in Eq. (66). Finally, for corners of the entanglement cut (where three surfaces intersect), the mapping is specified by

[TABLE]

where $i\neq j\neq k$. The above mappings are depicted schematically in Fig. 5.

It is straightforward to show that these mappings preserve all of the commutation relations between the terms in the EH defined by Eqs. (63) and (64). Considering only their support on the effective link variables, the hinge terms can be mapped onto the TFIM, with the exception that there are two (three on the corners) distinct $\mathcal{Z}\mathcal{Z}$-like Ising nearest-neighbour terms and, crucially, that the model is defined on the “cage” (wire-frame) topology instead of on a cycle. The link variables along the cage are coupled to the surface degrees of freedom by the fact that the $\mathcal{Z}\mathcal{Z}$-like terms along the cage anti-commute with the $\mathcal{X}_{i}$ surface terms. Recall, however, that the $\mathcal{X}_{i}$ terms in the surface TFIP model Eq. (67) occur at higher order in perturbation theory than the $\mathcal{B}$ terms, as indicated by the coeffecient $\epsilon$ of the transverse field term. Thus, to lowest order, we find that the hinge and corner contributions to the EH are mapped onto the TFIM acting on effective spin-$1/2$ degrees of freedom living on the edges of the cubic cut. In addition to distinguishing it from the $d=3$ toric code, the EH also clearly distinguishes the fracton phase from stacks of weakly coupled $d=2$ toric codes, for which the EH would consist of independent TFIMs along each row and column of the square lattice on every surface of the cubic cut. The distinction between the fracton phase and $d=3$ or decoupled stacks of $d=2$ topological orders is also evident in the “cage-like” nature of the EH. In the future, it would be interesting to study the evolution of the entanglement spectrum under the flux-string condensation discussed in Refs. Ma et al. (2017); Vijay (2017), whereby stacks of $d=2$ toric codes are strongly coupled to arrive at the X-cube.

Following this analysis, we can generalise to any entanglement cut given by a rectangular prism. The X-cube EH can be mapped onto the TFIP model along the surface of the cut, with the hinge and corner contributions mapping onto the TFIM along the cage formed by the edges of cut. Recall that for topologically ordered phases, we found that the EH is universally given by a $\mathbb{Z}_{2}$ invariant Hamiltonian acting on effective spin-$1/2$ variables for generic entanglement cuts. Given the geometric nature of fracton order, it should be no surprise that the EH is correspondingly sensitive to the geometry of the entanglement cut. To wit, the EH of the X-cube model can be mapped onto a TFIP model along the surface of the cut and to the TFIM along the hinges of the cut only for cuts which respect the $d=2$ planar subsystem symmetry of the X-cube. For instance, for a cut along the [111]-direction the EH will not take this simple form, and so generically, the EH will not fall into the class of subsystem symmetric invariant Hamiltonians acting on effective qubits—it is only for cuts respecting the $d=2$ planar subsystem symmetry that the EH will be invariant under $d=1$ subsystem symmetries.

Returning to a cubic entanglement cut, there is an ambiguity in the mapping of terms in the EH due to the self (or Kramers-Wannier) duality of the TFIP model, similar to what we observed for the $d=2$ toric code. This duality is related to our choice to place effective degrees of freedom on the dual lattice on the surface whereas we could equally have placed them on the vertices of the original lattice 888We choose to use the dual lattice description here since this makes the hinge and corner mappings appear more natural.. As with the $d=2$ toric code, the duality for the X-cube EH becomes exact once we consider the NLSS constraints, whereas for a two-dimensional TFIP model the duality is only exact up to the subsystem NLSS constraints along any rigid line. For the effective TFIP model onto which the EH maps, these constraints are primarily responsible for the fractonic behavior among the $\mathcal{B}_{p}$ terms in the EH (67). As discussed in Ref. Schmitz et al. (2018), all $6R+1$ independent topological NLSS for the cubic cut are guaranteed by the existence of a sub-extensive number of topological constraints. The constraints amongst the cubic stabilizers are generated by the product of all cubes which contain vertices in any given plane perpendicular to a coordinate direction. Similarly, the topological constraints among the vertex stabilizers is generated by the product of all vertex stabilizers which are entirely supported in any coordinate plane. Thus the product of all stabilizers in $A$ for any one of these planes forms a ribbon NLSS as depicted in Fig. 4. Note that the cube ribbon parallel to the surface is supported on the edge where the ribbon bends. We can characterise the resulting NLSS constraints as

[TABLE]

where $\text{ribbon}_{\perp}$ is any rigid string of vertices wrapping around the entanglement cut along coordinate directions (thus corresponding to a ribbon perpendicular to the surface) and $\text{ribbon}_{\parallel}$ is a rigid ribbon of plaquettes and edges wrapping around the cut along coordinate directions. After the mapping onto an effective TFIP model, the constraints among the terms of the EH are given by

[TABLE]

The first (plaquette) constraint is naturally enforced whereas the transverse field constraint is not and must be enforced by hand; since this term occurs at higher order in perturbation theory, however, we need only focus on the first constraint. This constraint is responsible for a $\mathbb{Z}_{2}$ subsystem charge conservation which is a key signature of fractonic behaviour Schmitz (2018); You et al. (2018); Song et al. (2018), i.e., for the X-cube model, charge is conserved globally but also along each plane. For the TFIM along the hinges of the cut, notice that the product of all parallel ribbon NLSS along any one direction forms a cage. It is this NLSS that constrains the hinge TFIM such that

[TABLE]

so that after the mapping we have

[TABLE]

The $\mathbb{Z}_{2}$ charge conservation is thus enforced for the entire cage formed by the hinges of the cut and not just for the boundary of a single surface. Finally, we note that any parallel ribbon NLSS which straddles the outer edge of a surface satisfies the definition of a superficial NLSS as all the cut cube stabilizers in this product are along the surface. This lends credence to our assertion that superficial NLSS signal the presence of subsystem symmetries where the superficial NLSS is formed. However, the existence of superficial NLSS is more pertinent to the discussion of Haah’s cubic code and SSPT models, so we postpone a detailed discussion of superficial NLSS until Section III.3 and IV.

III.3 Type-II Fracton Order: Haah’s Cubic Code

We now study the entanglement spectrum of a type-II fracton model, exemplified by Haah’s cubic code Haah (2011). The model is defined on a cubic lattice with two qubits living on each vertex of the lattice. The model is described by the Hamiltonian

[TABLE]

which consists of two stabilizer types, both of which are associated with every elementary cube of the lattice. The first type is composed of X-type Pauli operators, $G_{c}^{X}$, and the other is composed of Z-type Pauli operators $G_{c}^{Z}$, with their precise form shown in Fig. 6. The curious properties of this model are reviewed in Ref. Nandkishore and Hermele (2018).

We start by considering an $R\times R\times R$ cubic entanglement cut along the coordinate directions of the lattice, so as to contrast type-II fracton order with type-I. Given the form of the stabilizers as shown in Fig. 6, we require at least third-order perturbations before any contributions arise in the EH along the surface of the cut. These contributions result from the $(3,5)$ or $(5,3)$ cut stabilizers, coming either from $G_{c}^{X}$ or $G_{c}^{Z}$ depending on which surface we are considering. Regardless of which surface we consider, away from the hinges of the entanglement cut there are two types of terms in the EH for each plaquette of the square lattice on the surface. One such example is depicted in Fig, 7(a), with all others unitarily equivalent. Contributions to the EH from the surface of the cut $s$ are hence give by

[TABLE]

where $\tilde{B}_{p}^{X}=\left(G_{c}^{x}\right)_{A}$ or $\left(G_{c}^{x}\right)_{B}$ depending on the particular surface under consideration, and likewise for $B_{p}^{Z}$.

Along the hinges of the entanglement cut, we find contributions at first, second, and third order coming from the $(1,7),(2,6)$ and $(3,5)$ cut stabilizers which intersect along the hinge (see Fig. 7). The hinges contribute terms of the form:

[TABLE]

where $\tilde{B}_{e}^{X(Z)}$ are defined the same way as the surface terms. The precise form of the cut stabilizer depends on whether the $XX$, $ZZ$, or $II$ corner of the $G_{c}^{X},G_{c}^{Z}$ stabilizers are included in the cut edge. For a given edge, if an $II$ corner is included along the edge of one stabilizer type, the $XX$ or $ZZ$ corner is included along the edge of the other stabilizer type, as can be easily checked.

Unlike all cases considered heretofore, there does not appear to be a familiar Ising-like model onto which the surface EH maps. In fact, we could have anticipated that would be the case given the discussion regarding the role of subsystem symmetries in the our analysis of the X-cube model, where the surface EH maps onto the TFIP model only for entanglement cuts along the planar subsystem symmetries. Similarly, we expect that the EH for cuts respecting the fractal subsystem symmetry of Haah’s code will be mapped onto a subsystem symmetric Ising-like model—this is indeed the case for a [111] cut, as we discuss later in this section.

Along the hinges of a cubic entanglement cut however, we do find mappings of the EH onto familiar models as long as we ignore any commutation with the surface terms. By considering the mutual commutation relations between the hinge terms specified in Fig. 7(b), we find that these contributions to the EH can be mapped onto the $d=1$ TFIM. This is the case for the six hinges which contribute $(1,7)$ and $(3,5)$ cut stabilizer terms. In contrast, contributions from the remaining six hinges coming from $(2,6)$ cut stabilizers can be mapped onto a CSS version of the $d=1$ cluster model i.e., onto the terms given in Fig. 7(c) 999To see this, coarse-grain the $d=1$ chain such that there are two qubits per unit cell. After this coarse-graining, there remain only two stabilizers types which are no longer equivalent, up to translates. These can then be mapped via a local Clifford circuit consisting of an operator for each unit cell, such that the resulting stabilizer code corresponds to Fig. 7.. Finally for the cubic cut, we note in passing that along both the surface and the hinges, there exists a self-duality whereby all $X$-type terms can be exchanged with the $Z$-type terms, leaving the EH invariant.

Although all NLSS are not well-understood for Haah’s code, for the cubic cut there exist $12R-2$ independent NLSS Schmitz et al. (2018). Fourteen of these can be generated using the topological constraints found in Ref. Schmitz (2018) 101010The “star” pattern in Fig 5a of Ref. Schmitz (2018) has two independent versions given by exchanging the configurations in the three $[111]$ layers of the triangular lattice over which the pattern is periodic; see the reference for more details. For example, two are implied by the topological constraints given by the product over all stabilizers of a single type, which enforces

[TABLE]

Besides the remaining independent NLSS implied by the other twelve known constraints, one can also show that there are no superficial NLSS along any single surface of the cubic cut.

As anticipated, the cubic entanglement cut does not reveal much about the entanglement structure of Haah’s code since it does not respect the (fractal) subsystem symmetries, which played a significant role in our analysis of the X-cube model. Indeed, there exist more illuminating cuts for which we find contributions to the EH at first and second order along the surfaces of the cut. In particular, this is the case for flat entanglement surfaces perpendicular to the $[\pm 1,\pm 1,\pm 1]$ directions, as these cut along the corners of the cubic stabilizers. As a consequence of the three-fold rotational symmetry about the $[111]$ direction and previous results discussed in Ref. Schmitz (2018), the most interesting case is that of the surface perpendicular to the $[111]$ direction. Specifically, let us consider an $R\times R\times R$, $45\degree,45\degree,60\degree$ parallelepiped, where the plane of the $60\degree$ angle is perpendicular to $[111]$ and half of the parallelepiped unit cell forms a corner of a unit cube, as depicted in Fig. 8. As four of the surfaces for this cut correspond to faces of a cubic cut, the contribution to the EH from those faces is the same as that in Eq. (75). However, each of the $[111]$ surfaces form a triangular lattice and have contributions distinct from those discussed prior. We look at the surface farthest from the orgin, with results for the closer face obtained by swapping stabilizer types. This surface contains $(2,6)$ $Z$-type and $(6,2)$ $X$-type cut stabilizers, such that its contributions to the EH are given by:

[TABLE]

where $\widetilde{TT}_{\triangle}=(G_{v}^{X})_{A}$ is formed by two stacked triangle operators as shown in Fig. 8 and where $\triangle$ corresponds to the vertex $v$ just above the triangluar plaquette in the $[111]$ direction. Here, the sum is restricted to only one set of triangles, i.e., it only goes over upward-pointing triangles.

Within the boundary layer of $A$, the $TT_{\triangle}$ operator forms a single triangle operator, strongly reminiscent of the Newman-Moore model Newman and Moore (1999). This model is defined on a $d=2$ triangular lattice, with the Hamiltonian given by

[TABLE]

where the sum runs over all sets of nearest-neighbour spins $i,j,k$ living on the three vertices of one of the upward-pointing triangles. To make the connection of the $[111]$ surface EH to the Newman-moore model precise, we introduce one effective spin-$1/2$ degree of freedom for every two-qubit unit cell along the entanglement surface and map the terms as follows:

[TABLE]

where $\tilde{\mathcal{T}}_{\triangle}=\widetilde{{}\mathcal{Z}_{\triangle_{1}}\mathcal{Z}_{\triangle_{2}}\mathcal{Z}_{\triangle_{3}}}$ such that $\triangle_{i}$ is one of the effective spin-$1/2$ variables forming $\triangle$ and where $v=i$. The $[111]$ surface Hamiltonian is then mapped onto

[TABLE]

which we recognise as the Newman-Moore model Eq. (79) in the presence of a transverse field. Hence, in contrast with the cubic entanglement cut, for the cut depicted in Fig. 8 we find a natural mapping of the EH for Haah’s code onto a transverse field Newman-Moore model, acting on effective spin-$1/2$ degrees of freedom. Taken alongside our results for the X-cube model, this clearly illustrates the non-trivial geometric dependence of fracton phases, whose EH is not universally equivalent to some effective spin model, which is the case for topologically ordered phases. Instead, for both type-I and type-II fracton phases, it is only for specific entanglement cuts i.e., those compatible with the subsystem symmetries of the phase, that the EH maps onto a subsystem symmetric model. For Haah’s code, the mapping of the $[111]$ surface EH is onto the Newman-Moore model Eq. (81), which is invariant under fractal subsystem symmetries Yoshida (2013); Williamson (2016); Devakul et al. (2019). Our results hence illuminate the crucial role played by subsystem symmetries in the entanglement structure of fracton phases.

As before, we still need to consider the effect of the projection due to the NLSS. In Appendix C, we find that the number of independent NLSS for the cut considered in Fig. 8 are $12R-2$. Fourteen of these are related to the same topological constraints as discussed for the cubic cut and can be explicitly derived from them Schmitz (2018). Unlike a cubic entanglement cut however, for which there are no superficial NLSS along the surfaces, in this case there do exist superficial NLSS along the $[111]$ surface. Consider the $(3,5)$ cut $Z$-type stabilizer such that its support forms a triangle in the boundary layer of $A$. As with the Newman-Moore model, the product of triangle stabilizers forming a Sierpinski fractal only has support on the triangles forming the corners of the fractal. For a fractal of size of the $[111]$ surface, all support is removed from $A$, leaving such a stabilizer group member in $G_{B}$, as shown Fig. 9. This does not affect the terms in Eq. (78) since such terms contribute at third-order in perturbation theory and so can be safely ignored. Nonetheless, this does prevent any local first order contributions to the EH of the form $X_{v}^{1}I_{v}^{2}$ for a vertex $v$ in the triangular boundary layer of $B$. While such terms commute with all complete stabilizers in $B$, they do not commute with the superficial NLSS, which is a consequence of the fractal subsystem symmetries present in Haah’s code. As we discuss in detail in Sec. IV, although such local terms are disallowed due to the superficial NLSS, there exist quasi-local combinations of such operators which may yet survive, with consequences for the entanglement structure of the phase. We return to this point in Sections. IV and V.

As an aside, we note that forthcoming work by one of us Schmitz has found a coupled layer construction of Haah’s code, similar in spirit to the coupled layer construction of the X-cube model put forth in Refs. Ma et al. (2017); Vijay (2017). Much like the EH for the X-cube, which clearly distinguishes the fracton phase from decoupled stacks of $d=2$ topological orders, we expect that the EH for Haah’s code should likewise be distinct from decoupled stacks of the underlying $d=2$ SSPT layers. We leave an analysis of the effect of the coupling procedure on the entanglement structure to future work.

III.4 Edge-ES Correspondence

Following the original proposal by Li and Haldane Li and Haldane (2008), in which the low-energy entanglement spectrum was conjectured to match the low-energy spectrum of physical boundary modes, the edge-ES correspondence and variations thereof were found to hold in many gapped topological phases. For chiral topological orders in $d=2$, whose edge excitations are well known to be described by a 1+1D conformal field theory (CFT), it was shown that the reduced density matrix of a bipartition in the bulk is equivalent to the thermal density matrix of the CFT describing the low-energy boundary dynamics Qi et al. (2012); the edge-ES correspondence was similarly established for a series of fractional quantum Hall states in Refs. Regnault et al. (2009); Papić et al. (2011); Chandran et al. (2011); Dubail et al. (2012); Cano et al. (2015). A geometric proof later extended the applicability of the edge-ES correspondence to physical systems with an approximate Lorentz invariance at low energies Swingle and Senthil (2012), including topological insulators, which belong to the class of symmetry protected topological (SPT) phases. In SPT states, which are short-range entangled states with non-trivial protected boundary modes, the ES harbours a corresponding symmetry protected degeneracy Pollmann et al. (2010); Turner et al. (2010); Fidkowski (2010); Alba et al. (2012); Choo et al. (2018).

While both chiral topological orders and SPT states in $d=2$ host non-trivial gapless edge excitations, with corresponding signatures of these states appearing in the bulk ES, the boundaries of non-chiral topologically ordered states, examples of which include the $d=2,3$ toric code, are generically gapped. Nonetheless, it was shown in Refs. Ho et al. (2015, 2017); Koch-Janusz et al. (2017); Luo et al. (2018) that a version of the edge-ES correspondence holds for such phases as well. Specifically, Ref. Ho et al. (2015) studied the Wen-plaquette model, which is unitarly equivalent to the $d=2$ toric code, in the presence of arbitrary local perturbations. While they did not find an exact matching between the low-lying spectrum of the physical edge Hamiltonian and the low-lying ES for generic perturbations, as is true for chiral topological states, a “weak” edge-ES correspondence was found to hold. The weaker form of the edge-ES correspondence is encapsulated by the fact that both the EH and the edge Hamiltonian universally belong to the class of $\mathbb{Z}_{2}$ invariant $d=1$ Hamiltonians acting on effective spin-$1/2$ degrees of freedom. We now provide evidence that a similar edge-ES correspondence generically holds for the gapped fracton phases considered in this paper as well.

We first note that for a perturbed stabilizer Hamiltonian (6), the derivation of the edge Hamiltonian (the Hamiltonian describing a physical boundary of the system), proceeds analogously to that of the entanglement Hamiltonian i.e., through UPT. That is, once a specific physical boundary is specified, we can repeat the perturbative analysis developed in Sec. II.3 (or equivalently, in Ref. Ho et al. (2015)) in order to find $H_{\text{edge}}$ to any given order in perturbation theory. For all models considered in this section, we find that to lowest order,

[TABLE]

up to shifting and rescaling. Here, $\simeq$ implies that both models belong to the same class of Hamiltonians.

In agreement with the results of Ref. Ho et al. (2015), we find that both $H_{\text{edge}}$ and $H_{\text{ent}}$ for the $d=2$ toric code can be universally mapped onto a $\mathbb{Z}_{2}$ invariant, self-dual TFIM acting on effective spin-$1/2$ variables along the one-dimensional entanglement cut. Similarly, for the $d=3$ toric code, both the EH and the edge Hamiltonian can be generically mapped onto an effective $d=2$ TFIM which is invariant under a global $\mathbb{Z}_{2}$ symmetry. In fact, it known on general grounds that the boundary theory for a $\mathbb{Z}_{2}$ topologically ordered phase in $D$-dimensions is unitarily equivalent to a $\mathbb{Z}_{2}$ invariant Ising theory in $D-1$ dimensions Freed and Teleman (2018), where the $\mathbb{Z}_{2}$ symmetry is enforced by the bulk topological order. Along with our results in Sec. III.1, this clearly establishes the edge-ES correspondence for the toric code. Ref. Freed and Teleman (2018) additionally shows that a topologically ordered phase with electromagnetic duality in the bulk hosts a boundary theory which is Kramers-Wannier (or self) dual, which we demonstrated was also the case for the EH of the $d=2$ toric code.

For both type-I and type-II fracton orders, typified by the X-cube and the cubic code models respectively, we found that the EH cannot be universally mapped onto some Ising-like model. Instead, it is only for entanglement cuts along the subsystem symmetries (planar for the X-cube, fractal for Haah’s code) that the EH can be mapped onto an Ising-like model invariant under a sub-extensive set of $\mathbb{Z}_{2}$ subsystem symmetries. The geometric sensitivity of the entanglement Hamiltonian hence serves to clearly distinguish fracton order from topological order, a feature which carries over to the boundary Hamiltonian of fracton phases as well. Curiously, much as the bulk $\mathbb{Z}_{2}$ topological order enforces global $\mathbb{Z}_{2}$ invariance in the edge and entanglement Hamiltonians for the toric code, the bulk $\mathbb{Z}_{2}$ fracton order enforces a subsystem $\mathbb{Z}_{2}$ invariance on the edge and entanglement Hamiltonians for both the X-cube and the cubic code. For the EH, this invariance is enforced through the NLSS, while for the edge Hamiltonian we expect that it is a consequence of the Wilson/’t Hooft operators in the bulk. Based on our results, it seems reasonable that results established in Ref. Freed and Teleman (2018) for topological order can be extended to gapped phases with fracton order as well. In other words, we expect that the boundary Hamiltonian (equivalently, EH) for a system with fracton order can be mapped onto a subsystem symmetric Ising model, as long as the physical boundary (entanglement cut) lies along the planar (for type-I) or fractal (for type-II) subsystem symmetries present in the bulk. While boundary theories of fracton phases have yet to receive much attention, we note that they have been studied in some detail for the X-cube model in Ref. Bulmash and Iadecola (2018), whose results match ours where there is overlap.

IV Quasi-Local Entanglement in SSPT Phases

In this section, we look at the entanglement Hamiltonian of stabilizer codes describing subsystem symmetry protected topological phases You et al. (2018, 2018); Devakul et al. (2019, 2018a). In particular, we consider the $d=2$ cluster model, using a CSS variant of the coarse-grained version introduced in Ref. Williamson et al. (2018a)111111Our version differs from that discussed in Ref. Williamson et al. (2018a) by a local unitary which applies the Hadamard gate on the second qubit for every vertex.. This model is defined on the square lattice with two qubits assigned to each vertex and the Hamiltonian

[TABLE]

where $C_{p}^{X},C_{p}^{Z}$ are given by translates of the terms defined in Fig. 10.

In Ref. Williamson et al. (2018a) it was shown that this model, which is an example of a strong SSPT phase (a notion defined in Ref. Devakul et al. (2018a)), harbours contributions to the topological entanglement entropy for square entanglement cuts along the coordinate directions. Ref. Williamson et al. (2018a) dubbed such contributions “spurious” since they survive the usual subtraction scheme defining the TEE Kitaev and Preskill (2006); Levin and Wen (2006), which is expected to vanish identically in gapped phases without long-range entanglement i.e., without topological order. Unlike contributions to the TEE for LRE states, such as the toric code, which are independent of the entanglement cut, the spurious contributions for the cluster state instead only arise for square cuts. This suggests that these spurious contributions arise as a consequence of superficial NLSS associated with the subsystem symmetry protecting the state. These subsystem symmetries are given by the product of all stabilizers of a single type along any rigid line in either coordinate direction. The corresponding operators are

[TABLE]

where $(x,y)\in\mathbb{Z}_{L}^{2}$ indexes a vertex of the graph. We now focus on the EH for square cuts to understand the significance of the subsystem symmetries, which we argue endow the strong-SSPT phase with “quasi-local" entanglement.

Due to the lack of any rotational symmetry in this model, each straight edge of the square entanglement cut has slightly varying contributions to the EH. For the top and left boundary, we find $(3,2)$ cut $X$ stabilizers and $(2,3)$ cut $Z$ stabilizers, whereas for the bottom and right boundary we find $(3,2)$ cut $Z$ stabilizers and $(2,3)$ cut $X$ stabilizers. At the corners of the cut, there is one $X$ and one $Z$ cut stabilizer, which are either $(2,3)$ or $(1,4)$ depending on the corner. For each edge along the cut, we find that the EH maps onto the $d=1$ TFIM, which is also what we would find for a square entanglement cut for the $d=2$ toric code (see Sec. III.1). However, the transverse field terms of the effective TFIM do not get projected onto a single charge sector, unlike the analogous terms for the toric code EH.

Crucially though, the nearest-neighbour Ising terms of the effective EH satisfy a “zero-charge” condition. This is attributed to the superficial NLSS given by the product of cut stabilizers in the set defining one of the subsystem symmetries along that boundary. As a consequence of the subsystem symmetry, such NLSS are concomitant with any rigid entanglement boundary, with the subsystem invariance hence fundamentally responsible for the survival of spurious contributions to the TEE. Surprisingly, it appears that the subsystem symmetries are behaving much as topological constraints do, resulting in contributions to the entanglement entropy which are typically associated with LRE systems. This is in consonance with recent work suggesting that, as a consequence of subsystem symmetries being intermediary between global and local (gauge) symmetries, strong SSPT phases may have a robust ground state degeneracy on topologically non-trivial manifolds May-Mann and Hughes (2018). Hence, strong SSPT phases are distinct from their SRE SPT counterparts and also from LRE topologically ordered phases, instead harbouring a pattern of “quasi-local” entanglement (QLE) in their ground states.

More precisely, we can view the subsystem symmetries as enforcing a psuedo-constraint in the following way: suppose we impose a restriction on the kinds of Pauli operators which are allowed to act on our system. Viewed as a limitation on the types of perturbations which are allowed, these must form a subspace of the Pauli space such that products of any of the restricted operators also belong to the set of restricted operators. In this case, let us now consider all operators which commute with the subsystem symmetries of the $d=2$ cluster states (modulo any of the stabilizers). Such a subspace $\mathcal{R}\subseteq P$ is generated by the set $\{XI,IZ\}$ for every vertex of the square lattice. Due to the restriction, not all parts of the stabilizers matter, i.e., we can project out the parts of the operators which commute with everything in $\mathcal{R}$, which happens to be $\mathcal{R}$ itself. As a result, the stabilizer set becomes two commuting copies of the Ising plaquette model and, significantly, the subsystem symmetries are projected onto the identity after the restriction. This implies that each subsystem symmetry imposes a psuedo-$\mathbb{Z}_{2}$ conservation law once the restriction is enforced. This restriction is also reflected in the EH.

Nïavely, one might think we have left out first-order contributions to the EH coming from $IX$ and $ZI$ perturbations. But while these terms locally commute with all complete stabilizers, they do not commute with the subsystem symmetry and thus do not commute with the superfical NLSS. However, such terms can contribute if, upon inclusion of their coefficients, they collectively commute with the superficial NLSS i.e., act quasi-locally. In other words, we can construct an operator $\mathcal{O}=\sum_{i}\chi_{i}\mathcal{O}_{i}$ which commutes with the superficial NLSS despite being composed of local terms $\mathcal{O}_{i}$ which do not individually commute with the superficial NLSS. In order for this to occur, the coefficients $\chi_{i}$ must conspire in a manner which preserves the invariance of $\mathcal{O}$ under the subsystem symmetries, implying a measure of non-locality in the system. In the language of Ref. Williamson et al. (2018a), the existence of such an operator $\mathcal{O}$ implies the absence of a linearly-symmetric local unitary transformation which can transform the cluster state into a trivial product state, thereby demonstrating the strong SSPT nature of the cluster state. In fact, the quasi-local nature of the strong SSPT phase is implicit in the classification developed in Ref. Devakul et al. (2018a): the unitary circuit required to transform a strong SSPT into a trivial product state, while of finite depth, must have support over a sub-extensive number of sites i.e., is “infinite-length” in the thermodynamic limit. This is in contrast with both SPT states, for which the corresponding unitary would have support over a finite number of sites, and with topologically ordered states, for which it would need support over an extensive number of sites.

Thus, we see that strong SSPT states are quasi-locally entangled, since they possess features characteristic of long-range entangled states, but only for entanglement cuts and perturbations respecting the subsystem symmetry: for a square entanglement cut, as long as the perturbations respect the subsystem symmetry, the EH is equivalent to the $d=1$ TFIM. If we allow perturbations which break these “gauge-like” symmetries Nussinov and Ortiz (2009)121212A gauge-like symmetry is one which acts independently on a $d$-dimensional sub-region of a $D$-dimensional system, with $0<d<D$., we generically find additional non-local terms in the EH, destroying its mapping onto the TFIM. In contrast with the EH for the $d=2$ toric code, which is equivalent to the $d=1$ TFIM for generic perturbations, the EH for the $d=2$ cluster state only maintains this equivalence for a certain class of perturbations; equivalently, the mapping of the EH onto the $d=1$ TFIM is protected by the subsystem symmetry. Hence, we see that the EH for the cluster state resembles that for a topologically ordered state but only as long as the subsystem symmetry is preserved. Since the entanglement entropy can be extracted from the EH, this further implies that the so-called spurious contributions to the TEE in the strong SSPT phase are in-fact relaying the quasi-local entanglement present in the state. It is in this sense that we assert that these seemingly ersatz contributions to the TEE—although not topological in nature—possess a measure of robustness beyond that found in conventional SPT phases and are indicative of the gauge-like nature of subsystem symmetries.

V Discussion: Role of NLSS

The primary result of Section II, which is the derivation of the EH Eqs. II.3 and II.3 for stabilizer codes in the presence of arbitrary (weak) local perturbations, leads to several interesting entanglement Hamiltonians, as we have established in the preceding sections. As a corollary, we have also provided evidence that a weak edge-ES correspondence, known to hold for non-chiral topologically ordered states Ho et al. (2015, 2017); Koch-Janusz et al. (2017); Luo et al. (2018), also extends to systems with fracton order. Another key result, which we now expound upon, is the existence of non-trivial features in the entanglement structure which are protected by the NLSS.

Thus far, all of our results have been established for the ground states of the perturbed systems. To extend this, we conjecture that the EH for low-energy, perturbed eigenstates takes the same form as that for the ground state. More precisely, we posit that the EH for a low-energy eigenstate of the perturbed system is also a sum of projected cut stabilizers, but with coefficients which may differ from those for the ground state. The validity of this conjecture comes down to evaluating the validity of the UPT unitary acting on the low-energy states. That is, if we calculate Eq. (II.3) with $U$ acting on a low-energy stabilizer eigenstate $\ket{\mathbf{k}}$, the analysis proceeds identically (up to the exact coefficients) so long as the excitations do not lie near the boundary defined by the cut. An important consequence of this is that for the emergent Gauss law $\prod_{i\in F_{A}}O_{i}=\left(g_{\partial}\right)_{A}$ from Eq. (17), Eq. (48) becomes

[TABLE]

where $F_{\partial}$ is the set of stabilizer indices such that $g_{\partial}=\prod_{i\in F_{\partial}}O_{i}$ and $c$ represents the number of excitations contained in $F_{A}$ modulo 2 for the state $\ket{\mathbf{k}}$. This implies the surface charge of the EH is confined to the $c$-charge sector. If this holds, the charge of the surface EH can measure the charges in the bulk, which constitutes a “charge” bulk-edge correspondence. Here, the charge is defined by a global conservation law enforced by a topological constraint, which also enforces the existence of an NLSS (see Refs. Schmitz (2018); Schmitz et al. (2018)).

However, all excited states of the models considered here can have additional non-topological degeneracy i.e., a degeneracy which is not related to logical operators. As a result, it is not clear that we can ignore the choice of basis within the excited degenerate eigenspaces131313In UPT, this translates into the application of a constant (in $\lambda$) unitary which is block diagonal in the degenerate subspaces, but has no matrix elements between different eigenstates of the logical operators. Because of this logical operator selection rule, the ground space remains unaffected by the change of basis operator and we are hence justified in ignoring any such operators in our original analysis.. An intuitive argument for why one should be allowed to safely ignore non-topological degeneracies is as follows: in the original definition of the stabilizer Hamiltonian Eq. (6), we allowed for arbitrary coefficients $\{J_{s}\}$. Aside from modifying the final coefficients appearing in the EH, allowing arbitrary coefficients $\{J_{s}\}$ does not affect the results of Sections II.2 and II.3, but can remove any non-topological degeneracies from the excited states. Stated differently, simply by choosing a subset of the coefficients $\{J_{s}\}$ to be negative, we can make any valid stabilizers’ eigenstate the unperturbed ground state uniquely, up to the topological degeneracy.

Assuming this conjecture holds, the topological constraint is fundamentally responsible for enforcing the correspondence between the surface charge sector and the number of bulk $\mathbb{Z}_{2}$ charges, evidence that such a charge bulk-edge correspondence is endowed with a measure of topological invariance. For topologically ordered (LRE) gapped phases, this bulk-edge correspondence has a strict topological invariance, in that it is robust against any local perturbations and appears universally for any entanglement cut. Moreover, this conjecture also implies that flux conservation presents clear signatures in the ES for $d=3$ topological orders, such as the $d=3$ toric code model. If the entanglement boundary has non-trivial topology, then one can measure a surface flux analogous to Eq. (85), corresponding to the number of fluxes enclosed or encircled by the cut. This is a consequence of the NLSS being associated with flux and charge topological constraints in the bulk. On the other hand, for strong SSPT phases such as the $d=2$ cluster state, the bulk-edge correspondence is robust only for entanglement cuts and perturbations respecting the subsystem symmetry. For gapped systems with fracton order, there exists a mesh of both kinds of invariance (topological and subsystem), since although the total number of charge types is sensitive to the geometry, a subset of these are topological in nature i.e., exist universally for arbitrary local perturbations and for generic cuts.

Where SSPT phases and fracton systems differ from topologically ordered phases is in the existence of NLSS which are not invariant under all equivalent bulk transformations of $A$, but are only invariant under transformations within the intersection of $A$ with the subsystem symmetries. This is evidenced by the existence of superficial NLSS when the cut aligns with the rigid subsystems. For SSPTs these subsystem symmetries are not strictly topological, but rather quasi-local in nature (see previous section), while for fractonic systems their nature may be either topological (as is the case for the X-cube superficial NLSS) or quasi-local (as is the case for the fractal superficial NLSS in Haah’s cubic code). Though the quasi-local superficial NLSS can not be associated with a bulk-edge topological charge correspondence, such NLSS can nevertheless be attached to a “pseudo-constraint” provided by the subsystem symmetry and a restriction on the Pauli space. Here, a restriction refers to a limitation on the kind of Pauli operators which are allowed to act on a subspace, whereby the type of excitations are limited as well. Restrictions can be formalised in the language of linear gauge structures, developed in Ref. Schmitz (2018). In this way, a subsystem symmetry which is contained in $G$ implies a psuedo-charge if there exists a restriction such that the subsystem symmetry is elevated to a constraint. Thus, under the restriction, the subsystem symmetry implies a conservation law, as demonstrated explicitly for the SSPT example in Sec. IV. It was in this sense that we asserted that the seemingly spurious contributions to the entanglement, despite not being topological in origin, are a consequence of the subsystem symmetry and imply constraints analogous to those typically encountered in LRE states. The presence of such NLSS endows the strong SSPT with a measure of quasi-local entanglement beyond that encountered in conventional SPT states but lacking the full topological protection associated with topologically ordered states. Thus, the lexicon of gapped quantum phases must be expanded to include QLE states as intermediaries between those with either short- or long-range entanglement.

VI Conclusions and Outlook

In this paper, we have discussed the entanglement spectrum for the ground state of a general stabilizer Hamiltonian, first by considering the unperturbed Hamiltonian and then by analysing the effects of weak, local perturbations via UPT. This led us to find a version of the edge-ES correspondence between the entanglement Hamiltonian in the bulk and the edge Hamiltonian along the boundary of the system. More importantly, we found universal entanglement features which are preserved under such generic perturbations as enforced by NLSS symmetries of the stabilizer Hamiltonian. These features take the form of emergent surface charges/fluxes, which, via the NLSS and the emergent Gauss’ law it implies, are related to charges/fluxes in the bulk. We interpret this as a type of “charge” bulk-boundary correspondence. We applied our general results to several examples, including conventional topological order as encapsulated by the $d=2,3$ toric codes, type-I and type-II fracton order as encapsulated by the X-cube model and Haah’s cubic code, and strong SSPT phases as encapsulated by the $d=2$ cluster model.

Whereas the emergent surface charges/fluxes were found to possess a robust topological invariance for conventional topological order, both fracton and strong-SSPT phases demonstrate a non-trivial dependence on geometry, corresponding to their subsystem symmetries. Our results hence establish the ES as a clear diagnostic for not only fracton order, but also for strong SSPT phases, with the latter characterised by a qualitatively distinct pattern of quasi-local entanglement. This suggests that SSPTs, despite not having the robust topological invariance associated with LRE gapped states, exhibit entanglement features beyond those expected for familiar SPT phases, as reflected in the invariance of their ES under subsystem preserving perturbations.

We also conjectured that the surface charges/fluxes could be used to measure the bulk charges/fluxes for perturbed excited states so long as certain assumptions about UPT hold for these states. In principle, one should be able to place bounds on the range of low-energy excited states for which we expect the conjecture to hold; we leave a such a detailed analysis to future work. Further, in this paper we have restricted our attention to understanding the ES of fracton and SSPT phases described by stabilizer code Hamiltonians. It would be interesting to go beyond the stabilizer formalism and study the ES for fracton orders not described by commuting-projector Hamiltonians. While the edge-ES correspondence has been established more generally for $d=2$ non-chiral, non-Abelian topological orders with gapped boundaries Luo et al. (2018), it remains to be seen whether similar results hold for $d=3$ non-Abelian fracton models with gapped boundaries Vijay (2017); Prem et al. (2018); Song et al. (2018). As another extension, we note that although commuting-projector Hamiltonians cannot describe chiral phases Kitaev (2006), it has been suggested that gapped chiral fracton phases can instead be captured by the language of higher-rank tensor gauge theory Pretko (2017c); Prem et al. (2018). In the future, it would be interesting to study whether a stronger edge-ES correspondence, similar to that for chiral topological orders, exists for chiral fracton phases as well. More generally, understanding the structure of entanglement in gapless fracton phases is an important future direction, albeit one which lies beyond the techniques developed in this paper.

Acknowledgements.

The authors are especially grateful to Rahul Nandkishore for inspiring discussions and prior collaborations on fractons. We also acknowledge stimulating conversations with Danny Bulmash, Meng Cheng, Michael Hermele, Wen-Wei Ho, Tom Iadecola, Michael Pretko, Hao Song, Ruben Verresen, Dominic Williamson, and Han Yan. AS is supported by the Air Force Office of Scientific Research under award number FA9550-17-1-0183. SJ acknowledges support by the U.S. Department of Energy, Office of Science, Basic Energy Sciences (BES) under Award number DE-SC0014415. AP is funded by a PCTS Fellowship at Princeton University.

Appendix A Unitary Perturbation Theory

We discuss our method for perturbatively constructing a unitary operator $U(\lambda)$ which relates eigenvectors of an unperturbed Hamiltonian to the dressed eigenvectors of a perturbed Hamiltonian. Our derivation, which is a variation on the usual Schrieffer-Wolff perturbation theory Schrieffer and Wolff (1966); Bravyi et al. (2011), is designed to preserve the unitarity of $U(\lambda)$ to all finite orders in perturbation theory and is closely related to the Wegner-Wilson flow discussed in the context of many-body localization Imbrie (2016); Pekker et al. (2017).

Consider the case where one has an unperturbed Hermitian operator $H_{0}$ acting on members of a finite-dimensional Hilbert space $\mathcal{H}$ with eigenvalues $\{E_{n}\}_{n\in I}$ and a corresponding complete orthonormal set of eigenvectors $\{\ket{n}\}_{n\in I}$. Note that this set is not unique in the case of a degenerate spectrum. Now, suppose we add a Hermitian perturbation $V$, weighted by a control parameter $\lambda\in\mathbb{R}$, such that our new Hermitian operator becomes

[TABLE]

We now wish to approximate the new eigenvalues $\{E_{n}^{\prime}\}_{n\in I^{\prime}}$ and eigenbasis $\{\ket{n^{\prime}}\}_{n\in I^{\prime}}$ of the perturbed Hamiltonian at any finite order in $\lambda$. As discussed above, we can achieve this by finding a unitary $U(\lambda)$ which maps $\{\ket{n}\}_{n\in I}$ to $\{\ket{n^{\prime}}\}_{n\in I^{\prime}}$. However, we must define the mapping between the two bases, in order for which we make the following assumptions:

$\ket{n^{\prime}(\lambda)}$ is an eigenvector of $H(\lambda)$, 2. 2.

$\ket{n^{\prime}}$ is analytic within some non-zero radius of convergence about $\lambda=0$, 3. 3.

$\lim_{\lambda\to 0}\ket{n^{\prime}(\lambda)}=\ket{n^{\prime}(0)}=\ket{n}$, and 4. 4.

for any other $\ket{m}\in\{\ket{n}\}_{n\in I}$, its associated function $\ket{m^{\prime}(\lambda)}$ satisfies $\braket{n^{\prime}(\lambda)}{m^{\prime}(\lambda)}=\delta_{nm}$, for all $\lambda$ within the radius of convergence.

Taken together, these conditions imply that there exists an analytic, operator-valued function of $\lambda$, $U(\lambda)$, which transforms the unperturbed eigenbasis into the new eigenbasis. Further, this implies the existence of an analytic skew-Hermitian operator-valued function of $\lambda$, $F(\lambda)$, such that $U(\lambda)=\exp(F(\lambda))$. We now define the operator-valued function

[TABLE]

which is also analytic. We can expand both sides of this equation using the expansions

[TABLE]

and match the two sides, order by order. Rearranging Eq. (87), one has that

[TABLE]

where for any operators $A$ and $O$, $\text{Ad}_{A}O=[A,O]$. Now we can expand both sides in $\lambda$ and equate coefficients at each order in $\lambda$. The LHS is already expanded, whereas the RHS is a bit more complicated:

[TABLE]

where $\|\vec{a}\|=\sum_{i}a_{i}$ and where the super operator $\mathcal{F}^{(n)}$ is defined as acting as the identity for $n=0$ and, for $n>0$, having the action on any operator $O$:

[TABLE]

We then find the following relations, order by order:

[TABLE]

where $\mathcal{L}:=F^{(1)}$. To solve for each operator, we define the super-operator $\int_{H_{0}}$, with an action on an operator $O$ of

[TABLE]

The super-operator $\int_{H_{0}}$ is a partial inverse to $\text{Ad}_{H_{0}}$ up to an arbitrary operator which commutes with $H_{0}$. Furthermore, the kernel of this function is all operators commuting with $H_{0}$ which, by definition, includes all $D^{(n)}$. Hence, at order $n$, we apply the $\int_{H_{0}}$ super-operator to both sides of Eq. (92), where we note that the $m=0$ term is $D^{(n)}$, and so it drops out from the sum. Likewise, the $m=n$ term in the sum is $\mathcal{F}^{(n)}H_{0}$, which includes the term $-\text{Ad}_{H_{0}}F^{(n)}$, which maps to $-F^{(n)}$ under $\int_{H_{0}}$. Since all other terms only involve lower-order operators, one can recursively solve for all $F^{(n)}$ and $D^{(n)}$. At lowest order, we find that

[TABLE]

as quoted in Eq. (21) in the main text for the stabilizer Hamiltonian. Note that all $F^{(n)}$ are not unique as we are always allowed to add an operator which commutes with $H_{0}$ at each step of the recursion. Above, we are implicitly choosing this operator to be zero.

Appendix B Canonical Basis for the Pauli Group

In this Appendix, we describe the notion of a canonical basis for the Pauli group over $N$ qubits. As the Pauli group has been “Abelianized,” we can now treat it as a vector space over $\mathbb{Z}_{2}$, where addition is given by the product of Pauli operators and scalar multiplication is given by the power of a Pauli operator. However, it appears that we have lost the commutation rules for Pauli operators. We can recover this by introducing the symplectic form $\lambda:P\times P\to\mathbb{Z}_{2}$, whereby for $p_{1},p_{2}\in P$, $\lambda(p_{1},p_{2})=1$ if the two operators anti-commute, and is zero otherwise.

We start by considering the trivial example of a canonical basis as given by $\{X_{i},Z_{i}\}_{i<N}$. Clearly, this is a basis for $P$, but to use it along with $\lambda$ to expand any Pauli operator (in analogous fashion to an othonormal basis in an inner product space), one observes that this basis has the properties $\lambda(X_{i},X_{j})=\lambda(Z_{i},Z_{j})=0$ and $\lambda(X_{i},Z_{i})=\delta_{ij}$. Thus, one can expand any Pauli $p\in P$ as

[TABLE]

If we wish to have a similar expansion in some other basis, we demand a similar set of requirements. We define a canonical basis as a spanning set $\{p_{i}\}_{i<N}\cup\{\tilde{p}_{i}\}_{i<N}$ such that $\lambda(p_{i},p_{j})=\lambda(\tilde{p}_{i},\tilde{p}_{j})=0$ and $\lambda(p_{i},\tilde{p}_{j})=\delta_{ij}$. These conditions imply that we have an expansion similar to Eq. (95) for every Pauli operator within this basis. As the stabilizer group $G$ is mutually commuting, any basis for it suffices as (nearly) one-half of a canonical basis, but the other half, a canonical dual, is only unique up to members of $G$. This also suffices to show the dimension of any mutually commuting subspace of $P$ is maximally $N$.

Appendix C Recoverable Information for the $[111]$ cut in Haah’s code

Here, we derive the recoverable information for an $R\times R\times R$, $45\degree,45\degree,60\degree$ $[111]$ parallelepiped entanglement cut, as illustrated in Fig. 8. We use the definition

[TABLE]

where the last equality is a consequence of using pure states. To compute $s$, we use $s=s_{A}=|A|-d_{G}$ and the fact that for bounded simple cuts, all of $G_{A}$ is generated by stabilizers in $A$ Ma et al. (2018b). Note that we have also removed the minimization over all possible generating sets (see Eq. (14)) as there are no trivial constraints for Haah’s code. So as to visually facilitate the derivation, we note that we can visualise this cut as a stack of triangular lattice layers, each of which is off-set with respect to the others such that if we start at a vertex and proceed upwards, we pass through a triangle of the layer immediately above and an oppositely oriented triangle in the layer above that one. This combination of a vertex and two triangles is also one of the more illuminating ways of visualising the stabilizers of Haah’s code. If we consider the $X$ stabilizer as “upwards-pointing” since it follows the order: vertex $\leftarrow$ triangle $\leftarrow$ triangle, then the $Z$ type stabilizer is “downwards-pointing” since it follows the opposite order: triangle $\leftarrow$ triangle $\leftarrow$ vertex. Note that for each triangle forming a stabilizer, the top (bottom) triangular layer is only supported in one of the two qubits at a vertex, the middle triangular layer is supported on the remaining qubits, and the vertex layer is supported on both. We refer to the position of the stabilizer as that of the single vertex.

Clearly, the number of qubits in $A$ is $|A|=2R^{3}$. As for complete stabilizers, the top two layers contain no complete X stabilizers while the bottom two layers contain no Z complete stabilizers. In every other layer, all stabilizers are complete except for along the two edges near an acute angle of the parallelepiped. This implies $d_{G_{A}}=2(R-2)(R-1)^{2}$ such that $s=8R^{2}-10R+4$. As for cut stabilizers, there are clearly $8R^{2}-10R+4$ coming from all vertices in $A$ which were not counted while calculating $d_{G_{A}}$. Every triangle in the bottom layer, as well as all those along the edge, represents a cut $X$ stabilizer located in the two layers just below the surface and likewise for the top layer and cut $Z$ stabilizers. This contributes $4R(R+2)$ cut stabilizers. For the remaining layers, every upwards-pointing triangle along an obtuse-angle boundary of the parallelepiped is an additional cut stabilizer. This contributes $2(R-1)(2R-1)$. Applying the definition, we find that

[TABLE]

Curiously, this is the same value as that for a cubic cut.