Entanglement Entropy, Quantum Fluctuations, and Thermal Entropy in Topological Phases

Entanglement entropy in topologically ordered matter phases has been computed extensively using various methods. In this paper, we study the entanglement entropy of topological phases in two-spaces from a new perspective---the perspective of quasiparticle fluctuations. In this picture, the entanglement spectrum of a topologically ordered system is identified with the spectrum of quasiparticle fluctuations of the system, and the entanglement entropy measures the maximal quasiparticle fluctuations on the EB. As a consequence, entanglement entropy corresponds to the thermal entropy of the quasiparticles at infinite temperature on the entanglement boundary. We corroborates our results with explicit computation in the quantum double model with/without boundaries. We then systematically construct the reduced density matrices of the quantum double model on generic 2-surfaces with boundaries.


Introduction
Matter phases with intrinsic topological orders, or topological orders for short, are exotic, gapped phases of matter beyond the Landau-Ginzburg paradigm [1][2][3][4][5][6][7]. They not only have expanded our knowledge of phases of matter but also have important applications, including robust quantum memories and topological quantum computers [5,8].
Topological orders exhibit long range entanglement [9,10], which can be captured by the topological entanglement entropy (TEE) [7,11] of the system. Suppose a (2 + 1)dimensional topologically ordered system is virtually divided into two regions A and B, the entanglement entropy (EE) of the system is the sum of an leading area term, proportional to the length of the boundary between the two regions, and a universal subleading term proportional to log D, where D is the total quantum dimension 1 . As such, TEE can be used to distinguish topological orders unless certain topological orders has the same total quantum dimension.
There has been extensive works computing the EE in topological orders [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25]. In this paper, however, we propose a novel perspective of the EE in topological orders in two spatial dimensions, namely, the picture of quasiparticle fluctuations. Consider a bipartite system in a ground-state, e.g., as in Fig. 1(a). We show as one of our main results that the reduced density matrix ρ A is identified with the quasiparticle fluctuations on the entanglement boundary (EB) separating A and B, and the corresponding EE explicitly counts the allowed states describing quasiparticle fluctuations on the EB. We shall call such states EB-states. An EB-state is allowed when it ensures that no quasiparticles exist anywhere in the system except on the EB and that the total charge of the system and that of the EB are trivial. Hence, these quasiparticles cannot be seen by any global observer, say Gabe, of the system, as each of them is confined in a pair with its anti-quasiparticle on the other side of the EB. Such quasiparticle-antiquasiparticle pairs, which happen to be cut by the EB, comprise the quasiparticle fluctuations along the EB. We reckon that our notion of quasiparticle fluctuations complies with the usual notion of vacuum fluctuations in quantum field theory, which refers to that the vacuum consists of instantaneous particleantiparticle pairs everywhere. Particles in vacuum fluctuations cannot be observed unless the particle-antiparticle pairs are torn apart locally by sufficient amount of energy.
Entanglement entropy counts the EB-states because of the following. To compute the EE, Gabe has to trace out region B. Tracing out region B would blind Alice from knowing the existence of region B at all. What Alice would see is a world like that in Fig.  1(b), complete darkness (no quasiparticles in the bulk) with a glowing physical boundary (quasiparticles residing on the boundary). Alice will detect quasiparticles at the EB, which now appears to her as a physical boundary.
Furthermore, Alice is able to measure the EB-states on the EB and find that they are all equally probable because the reduced density matrix ρ A can be made a projector. She can then compute the thermal entropy (TE) of the boundary quasiparticle system at temperature T = ∞. This TE would coincide with the EE computed by Gabe. The infinite temperature can be thought as what offers the energy to tear apart the particle-antiparticle pairs at the EB.
To make the above picture of quantum fluctuations concrete and precise, we consider in particular the quantum double (QD) model of topological orders in two-spaces and focus on the ground-state Hilbert space of the model. The QD model [7,26,27] is a lattice 1 For a topologically ordered system with n type of anyons and di is the quantum dimension of the type-i Hamiltonian extension of the Dijkgraaf-Witten gauge theory [28] with discrete gauge groups. We focus on ground-states because anyon excitations in topological orders add no more physical content to the EE. In the QD model, we construct two useful types of bases-the holonomy bases and the fusion bases. The holonomy bases and fusion bases not only simplify the EE computation but also manifest the picture of quasiparticle fluctuations of EE. We compute the EE on a sphere in both bases, which leads to the same result, which agrees with the known result. This result corroborates our picture of quasiparticle fluctuations. We then extend our study to topological orders with gapped physical boundaries (PBs), separating the system from the vacuum, as opposed to the EB separating two regions of the system. See Fig. 2. We consider the extended QD model [29,30], dedicated to the scenario with PBs. With minor modifications, the holonomy and fusion bases work in this scenario, and hence the picture of quasiparticle fluctuations still applies. To be specific, we considered two representative cases, namely, a cylinder with region A being a meridian strip and a cylinder with region A being a longitudinal strip. 2 Using the holonomy and fusion bases, we offer closed-form formulae for the EE in both cases and an explicit example when the bulk gauge group is the permutation group S 3 . Such a formula in terms of the input data of the QD model has not been obtained before. We then systematically construct the reduced density matrices for topological orders on generic two-surfaces with gapped PBs.
We note that the EE on topological orders with PBs have also been computed recently on the extended QD model [22] and by using conformal field theory techniques [31,32].
We reckon that the fusion bases we construct are related to the fusion bases used in a recent work studying the generalized anyon exclusion statistics [33,34]. Hence, we expect to relate the EE in topological orders to the generalized anyon exclusion statistics. This expectation is reasonable because exclusion statistics is bound to the notion of thermal entropy. We shall leave this connection for future studies. This paper is organized as follows. Section 2 reviews the QD model and develops the holonomy bases. Section 3 constructs the fusion bases and explains the picture of quasiparticle fluctuations. Section 4 studies the cases with PBs, with a concrete example for G = S 3 . Section 5 discusses certain subtleties. The appendices collect certain calculations too detailed to be part of the main text.

Quantum double model and holonomy bases
In this section, we shall discuss the minimal and complete set of observables in a lattice topological gauge field theory. We will also discuss the relation between the gauge fields and these observables. This relation will enable us to study the entanglement from the perspective of observables.
We define a discrete gauge field theory on a 2D directed graph, with a finite gauge group G. The degrees of freedom (d.o.f.) are group elements on the edges. The Hilbert space is spanned by all possible configurations of these group elements on the edges. We identify the two states, (2.1) where on the RHS the group element is inversed and the arrow on the edge is reversed. For simplicity and illustration, we show only four-valent vertices but in general the spatial graph could contain any multivalent vertex. Then we define a gauge transformation at a vertex v by which is a left multiplication of h on the edges coming into the vertex v. A lattice gauge theory is defined by a Hamiltonian that commutes with any gauge transformation on the lattice. The quantum double (QD) model, also known as the Kitaev model [7], is such a lattice gauge theory, whose Hamiltonian reads where A v is a projection operator where δ abcd = 1 if abcd = 1, the unit element of G, and δ abcd = 0 otherwise. Here abcd is the holonomy around the plaquette p. Hence, B p projects onto states with trivial holonomy at p. In particular, when G = Z 2 (known as the toric code model), the Hilbert space is spanned by spins on the edges, and the Hamiltonian terms are where e represent the edges. Indeed, compared to Eq. (2.2), e into v σ x e is (the nonidentity) Z 2 gauge transformation with the Z 2 represented by {1, σ x }. The delta function in Eq. (2.5) becomes 1 2 1 + e around p σ z e .

Observables
Since all the operators A v and B p commute with each other, they form commuting observables of the toric code model, and their common eigenvectors form the basis states of the Hamiltonian. As a results, there is a Z 2 charge excitations at a vertex v in a state |ψ if A v |ψ = − |ψ , and a Z 2 flux excitation at a plaquette p if B p |ψ = − |ψ . The name "charge" and "flux" can be justified by the Aharonov-Bohm effect as follows. For a given edge e, the σ z e serves as the hopping operator that moves the charge from one end v 1 of e to the other end v 2 because A v 1,2 σ z e |ψ = −σ z e A v 1,2 |ψ . Hence, being the tensor product of σ z around the plaquette p, B p moves a charge around p by one turn and measures the Z 2 flux in the plaquette p. The names "charge" and "flux" thus follow.
On a torus, The minimal and complete set of commuting observables (also called the set of stabilizers) of the toric code model is given by where V is the total number of vertices, and P the total number of plaquettes. The two operators A V and B P are excluded from the above set because of the the two global constraints The operators Z 1 and Z 2 are the products of σ z along two arbitrary, non-contractible loops, one meridian and one longitudinal, of the torus. These two operators clearly measure the Z 2 fluxes going through the two holes of the torus (as seen in 3D space); they are responsible for the corresponding Aharonov-Bohm phases. We arrive at two bases of the Hilbert space. One is spanned by the spins on all the edges, while the other is given by the simultaneous eigenvectors of the operators in Eq.
(2.7). The latter corresponds to three types of physical quantities in a electromagnetic theory: Z 2 electric field strengths, magnetic field strengths, and Aharanov-Bohm phases due to nontrivial spatial topology. This establishes a duality between the observables and the gauge fields. We call the Aharonov-Bohm phases the topological d.o.f..

Holonomy bases: a non-local transformation
The above duality between the gauge fields and observables can be also made precise in the QD model of a generic group G. For our purposes in this paper, we focus on the states with trivial holonomy, i.e., B p = 1, everywhere. That is, we restrict ourselves to the subspace H Bp=1 of the total Hilbert space of the QD model.
We define a non-local transformation as follows. First, we pick an arbitrary vertex v 0 as the base point and choose a path on the graph from v 0 to v, for any other vertex v. Second, Given a configuration of the group elements a e on all the edges satisfying the trivial holonomy condition, we assign to each vertex v = v 0 a new group element where p vv 0 is the path chosen from v to v 0 . If the space has a non-contractible loop C α , we also assign a g α -the product of the group elements along a C α to C α . This way, we get a non-local transformation: Such a transformation is illustrated in Fig. 3. . The explicit transformation consists of g 1 = a 1 , g 2 = a 2 a 1 , g 3 = a 3 ,g 4 = a 5 a 3 , and g 5 = a 7 a 5 a 3 . Such a transformation from the a's to g's is clearly invertible.
For the quantum double model on a torus, g α are the holonomies along the two noncontractible cycles. See illustration in Fig. 4. In particular in toric code model case, they corresponds to the operators Z 1 and Z 2 .
The non-local transformation above defines a new basis of H Bp=1 , {|g v ; g α }, which we call a holonomy basis, as inspired by Ref. [35]. In a holonomy basis, an operator A v =v 0 acts as Figure 4. A basis of Hilbert space on torus, labeled by eigenvalues of A v , (v = 1, 2, . . . , V − 1), and Z 1 , Z 2 defined along the two non-contractible loops C 1 , C 2 .. There is a base vertex shown as dot in the diagram. The gauge charge at vertex each vertex v is formulated by the product of group elements along a path between the base point and v.
and A v 0 acts as (2.12)

Fusion basis and the quasiparticle picture of entanglement
In this section, we propose a new picture of the entanglement between the two parts of a bipartite topologically-ordered system. In this picture, we can compute the entanglement entropy by counting the configurations of quasiparticles on the boundary between the two parts that divide the system. We thus christen this new picture the quasiparticle picture. We shall continue our study using the QD model. Given a bipartite system, a state is entangled if it is not a tensor product of the states of the two subsystems. Entanglement entropy is a quantity that measures how much the two subsystems are entangled.
For a given state ψ in a bipartite system, the entanglement entropy is defined as where ρ A = tr B |ψ ψ| is the reduced density matrix. As such, the EE also measures one's ignorance of subsystem B, whose d.o.f. are traced out. One may also choose to trace out subsystem A but S E is independence of such choices. In this paper, we focus on the EE in the ground states of the QD model.

Entanglement entropy on a sphere in a holonomy basis
To motivate the quasiparticle picture and understand how the EE is computed in the new picture, let us reexamine a well-studied case, namely, the EE of the QD model with a finite group G on a sphere. In this case, the ground state is unique up to local unitary transformations, with trivial charge and zero flux everywhere on the sphere. This makes a holonomy basis handy for our computation of the EE. We virtually cut the sphere along an EB into a bipartite system. We would like to express the ground-state condition explicitly in terms of the charge d.o.f. on the EB. Doing this will manifest the contribution of quasiparticle fluctuations near the partition boundary to the EE and enable us to write down the reduced density matrix.
As illustrated in Fig. 5, we choose a base point v A (v B ) in subsystem A (B). We assume that the EB consists of L vertices. Since as far as the ground-state is concerned we can restrict the system to the subspace in which B p = 1 everywhere, plus the sphere has no non-contractible cycles, the non-local transformation (2.10) would trade the original d.o.f. on the lattice edges for the new d.o.f. on the paths linking the base points to the vertices on the EB. There should be two base points, one in subsystem A and one in B. Otherwise, if there is just one base point chosen, say, in A, then there must be a path between any two neighbouring vertices on the EB, rendering all paths being closed. But B p = 1 everywhere, all such closed paths must be trivial, making no sense.
Consequently, for each vertex on the EB, we have a charge d.o.f. g v in A and g ′ v in B. A ground state should be gauge invariant at each v, and has zero holonomy everywhere. The latter condition reads In terms of the new set of variables, the Hilbert space is spanned by the basis states To obtain the ground state on the sphere, we can pick any basis state and act on it by all the vertex operators A v =v 0 and A v 0 , which will project the basis state to a ground state. That is, pick a generic |g 1 , . . . , g L ; g ′ 1 , . . . , g ′ L , then by Eqs. (2.11) and (2.12), we have The second equality above is due to the trivial holonomy condition after the action of the vertex operators, namely, where the second equality is due to the trivial holonomy condition (3.2) before acting the vertex operators. By renamingg i back to g i again, we have a ground state where the second row follows the first row by an obvious suppression of indices. One can check that |Φ is indeed a ground state. Since We check that ρ A is a projection operator The normalization factor 1 |G| 2 in Eq. (3.6) is chosen such that ρ A is exactly a projector. Note that since ρ A must represent a mixed state, if ρ A can be made a projector, it means that all probable pure-states are equally weighted in ρ A . This leads to S E = −tr(ρ A log ρ A ) = log trρ A . Now that 8) and the EE computes as This result agrees with the known result (See for example Ref. [22]). In the above we defined a quantity W := trρ A . Clearly, Since all probable pure-states are equally weighted, W is simply the number of probable pure-states. The questions are, what are these pure states and how we may interpret them? To answer the questions, it would be better to Fourier transform the holonomy bases into what we call the fusion bases.

Fusion bases
We have just derived the reduced density matrix and computed the corresponding EE of the QD model on a sphere. This is done in a holonomy basis proposed in Section 2.2.
In the holonomy basis, the d.o.f. in either region A or B all live on the EB separating the two regions, and the d.o.f. in A and those in B group into entangled pairs, as seen in Eq. (3.5). This observation complies with the common sense that the entanglement between two regions actually occurs only near the EB between the two regions. It is also often said that the entanglement entropy measures the ignorance of an observer in region A (B) about region B (A). The computation in the holonomy basis not only manifests this common sense but also infers a new picture of entanglement, in which the ignorance can be made concrete in terms of the quasiparticle fluctuations on the EB. Imagine an observer in A. Tracing out region B means that the observer in A is completely refrained from knowing the existence of B, instead, the observer senses a boundary of A, which is the EB that now appears physical to the observer in A. Tracing out B must leave certain residue on the EB. Recall that the entire system is in a ground state and that the total state (3.5), which is pure, is Schmidt decomposed in the holonomy basis as superposed entangled pairs. Since now the observer in A cannot see region B at all, he would likely see quasiparticle excitations on the EB. These excitations are paired up with those as appeared to an observer in region B blinded from seeing region A, such that the entire system remains in a ground state. In the QD model, however, (G-charge type) quasiparticle excitations are better expressed in a basis not in terms of the group elements but in terms of the representations of the gauge group. As to be seen shortly, such a basis is a Fourier transform of the holonomy basis. Hence, what the observer in A sees would not be a single type of quasiparticles but in fact a mixed state of all possible quasiparticle states on the EB. In other words, the observer sees quasiparticles fluctuating on the EB. In the following, we will apply the non-local transformation (2.10) and rewrite the reduced density matrix in terms of quasiparticle d.o.f..
We shall Fourier transform the holonomy basis in terms of group elements to one in terms of the irreducible representations of the gauge group. We label the charges by the irreducible representations. This will result in a more transparent picture of the quasiparticle fluctuations near the EB. Define where j labels the unitary irreducible representations of G, and mn are matrix indices of the representations. The inverse transformation reads The above is the transformation on a state labeled by a single group element. The two transformations are unitary and inverse to each other, as can be proved by Peter-Weyl theorem in group representation theory. A state on the right hand side of the above equation can be regarded as a quasiparticle carrying the charge labeled by the irreducible representation. In the state (3.5), however, the holonomy basis consists of a chain of group elements for each of the two regions. Each such chain which would be transformed as in Eq. (3.11), resulting in a tensor product of states labeled by irreducible representations.
The tensor product can be interpreted as the fusion of a chain of charge quasiparticles, leading naturally to a fusion basis |j v m v ; η of the states, depicted as follows. Note that the fusion basis (3.12) shows that the quasiparticles on the EB must all fuse to be a trivial charge. That is, the total quasiparticle charge of region A including the EB is trivial. In fact, the reduced density matrix (3.6) is an identity matrix on all quasiparticle d.o.f. on the EB, combined with the projection P G , where P G is the average of right multiplication on all g v by h. It enforces a G-symmetry. This global symmetry is a consequence of gauge invariance condition at the base point in A, which requires the total charge of all quasiparticles in A to be trivial. The global G-symmetry can be viewed as the broken gauge symmetry due to the appearance of the EB. As a result, the nonzero eigenspace of ρ A are the space of all states of quasiparticles on EB that are invariant under this global G-symmetry. The η is nothing but a set of good quantum numbers of this global symmetry.
In the fusion basis, Eq. 3.6 becomes where η, η ′ label the fusion basis with j v m v fixed. Since the Fourier transform cannot prevent ρ A being a projector, the EE remains where W is the dimension of the space spanned by the fusion basis. Directly counting the number of the basis vectors yields (3.16) 3 We note that such a basis was adopted in understanding the entanglement entropy in lattice gauge theories by one of us in Ref. [19].
where use of j N jk l d j = d k d l is made in the second equality, and j 1 (j L ) is renamed to η 1 (η L−1 ). See Appendix A for more details. This result agrees with the previous result (3.8) obtained in the holonomy basis. This alternative derivation however has a deeper meaning, as we now elaborate on.

Entanglement measures quasiparticle fluctuations
The fusion-basis derivation in the previous subsection manifests a picture of entanglement in view of quasiparticle fluctuations. In this picture, we are able to answer the questions raised in the end of Section 3.1.
In the language of gauge field theory, the ground state(s) of QD model satisfy two conditions-zero charge everywhere and zero flux everywhere. These two conditions determine the nonzero eigenspace of ρ A , which can be spanned by the holonomy basis or the fusion basis developed in the previous subsections. In the following, we analyze the consequence of these two conditions on the charge d.o.f. g v and g ′ v on the EB . Zero-flux condition requires the charge d.o.f. on both sides of the EB to be paired up to a global symmetry, as concluded in Eq. (3.2). Here, we can imagine that the EB is fattened a bit to be a narrow strip. In Fig. 5, this requires g v = g ′ v h for some constant h at all vertices v on the EB. We shall consider the global symmetry on these charge d.o.f. later and let us ignore the constant h at this moment. Locally, charges appear in pairs (and the pair goes across the fattened EB) at each vertex on the EB. We say that the charge quasiparticles in an EB-state satisfy the: boundary matching condition : at all vertices v on the EB. Zero-charge condition has two consequences. First, at each vertex v on the EB, charge pairs appear with equal weight. That is, locally, the ground state must read as g |g A ⊗ |g B because it is the only state invariant under gauge transformations at v. Hence, at each v on the EB, the ground state locally reads as equal-weight charge pair : (3.18) The equal-weight charge pairs leads to a picture of quasiparticle fluctuations. A quasiparticle fluctuation is a charge pair across the EB. By the equal-weight condition (3.18), in the fusion basis these quasiparticles j v m v (as well as their fusion multiplicity d.o.f. η) are paired across the EB. When the system is partitioned into two regions by the EB, however, local observers in a A (in a mixed state ρ A ) can detect quasiparticles with equal probability on the EB. Our picture coincides with the quantum fluctuation in a vacuum state of a quantum field theory, which refers to a particle-antiparticle pair locally excited by an amount of energy. In our case, no actual energy is injected into the system, and a global observer (in a pure state |Φ ) cannot detect any quasiparticles anywhere because all quasiparticles on one side of the EB are confined with their counterparts on the other side. Second, the gauge-invariance condition at the base point of region A imposes a global symmetry, implemented by the P G in Eq. (3.13). All states |g v A in A ought to be invariant under the global right-multiplications over G. This global symmetry can be viewed as the broken gauge symmetry in A. Since as appeared to a local observer in A the EB is a PB, this gauge symmetry breaking is understood as due to the charge condensation on the EB [36,37].
To sum up, the ground state(s) of the QD model meet three conditions on the EB: the boundary matching condition, equal-weight local charge pairs, and a global symmetry P G .
The advantage of the holonomy bases is that the d.o.f. are well allocated into two regions. By "well" we mean that all the three conditions on the ground state(s) are inherited by the EB states of region A (B). In the following, we reexamine the three conditions from the eyes of the observers within A.
In the sphere case, the reduced density matrix ρ A is a projector, implying that the considered ground state |Φ admits a Schmidt decomposition:  (3.17), together with the equal weight condition (3.18) on chair pairs, implies that all possible charge quasiparticles will appear with equal probability on the EB in A. In the sphere case, this leads to that each of the W basis-states |Ψ i A is equally probable.
The EB-states |Ψ i A form a basis of all allowed states (in A) of the charge quasiparticles on EB. We summarize EB conditions on the EB states of region A (B) as follows.
1. They ensure that region A (B) has no quasiparticles except on the EB.
2. Locally, charge quasiparticles appear with equal probability everywhere on the EB.
3. If there are non-contractible loops in region A (B), the holonomy along these loops is unit element of G.
4. The total quasiparticle charge of either the entire region A or B is trivial, i.e., the quasiparticles in each allowed EB-state must fuse to be a trivial charge in either A or B.
Condition 3 is due to the special choice of the ground state which has trivial holonomy along these non-contractible loops. In general, if we choose other ground states |Φ ′ , we should modify this condition such that any observableÔ along these non-contractible loops in region A (B) should evaluates to the same value in |Ψ i A (|Ψ i B ) and in |Φ ′ . Condition 4 may look somehow odd at first glance that we require a global symmetry in A and one in B at the same time. We need both since they are the broken gauge symmetries in A and in B. In the sphere case, the two global symmetries coincide with P G but in general they may differ as the spatial topology of the system becomes more complicated. We shall see such example in the next section and will discuss about this condition in more details in Sec. 4.3.
Eq. (3.19) reads as a boundary matching condition upon gluing regions A and B along the EB, such that the total state |Ψ is a ground state that is free of quasiparticles everywhere.
Therefore, the EB-states of A and those of B are what are being entangled, as they are paired up. They are maximally entangled because all the EB-states are equally probable in the Schmidt decomposition. This physical picture is illustrated in Fig. 6. Each EB-state specifies a configuration of charge quasiparticles on the EB. The EB-states of A and those of B must match, such that the total state of the entire system is a ground state.
In terms of quasiparticles, the entanglement entropy entropy can be interpreted as a measurement of quasiparticle fluctuations. Quasiparticles on one side of the fattened EB are paired with their counterparts on the other side, and "paired" means "entangled". In fact, such pairing is the only contribution to the entanglement entropy, as we can see in the formula S E = log W , which counts how many states of such quasiparticles are allowed on the EB. The fact that charge quasiparticles appear with equal probability on the EB is a consequence of a equilibrium state in region A with maximal quasiparticle fluctuations across the EB. We say an equilibrium state has "maximal quasiparticle fluctuations" if the energy threshold to excite any charge pairs on the EB becomes zero. (A ground state is an eigenvector of A v = 1 and reads locally as a superposition of equal-weight charge pairs as in Eq. (3.18). The partial trace tr B A v becomes an identity operator, implying that quasiparticles on the EB can be excited without energy penalty as seen by the observers in A.) Hence we conclude:

Entangled spectrum corresponds to all possible states with quasiparticle fluctuations on the EB, and Entanglement entropy measures maximal quasiparticle fluctuations on the EB.
In other words, entanglement entropy S E = log W counts the number of entangled EBstates of regions A and B. To compute S E amounts to the problem of state counting W of the EB states satisfying the EB conditions. The EE (3.9) consists of an leading area term L log |G| and a universal subleading term − log |G|. The leading term results from the fact that there are L vertices on the EB where quasiparticle pairs are allowed, and that there are |G| types of quasiparticle pairs at each vertex (in either basis |g or |jmn ). The universal subleading term − log |G| results from the global symmetry imposed by P G in Eq. (3.13), which says the quasiparticle fluctuations across the EB should preserve this global symmetry. The result agrees with our previous observation that the nonzero eigenspace of ρ A = P G is the space of all states of quasiparticles on the EB that are invariant under the global G symmetry.

2D entanglement entropy = 1+1D thermal entropy
The physical picture described above also demonstrates a correspondence between the EE and thermal entropy (TE). Suppose the total system is in a ground-state |Φ . A local observer restricted in A sees no existence of region B but can measure the charge quasiparticles on the EB. The observer finds that there are W equally probable quasiparticle configurations on the EB. Under the micro-canonical ensemble assumption, he computes the TE of the quasiparticle system on the EB as log W . This result coincides with the EE we, as a global observer of the total system, obtained earlier. The local observer in region A certainly does not know that the TE was a consequence of a hidden region B but thought it was due to the quasiparticle fluctuation on the EB (a physical boundary as appeared to the observer. Note that the TE should be understood as one at infinite temperature because at T = ∞, all the quasiparticle configurations contribute to the partition function equally. We believe that the above correspondence in the sphere case can be generalized to a generic surface with non-trivial topology. If the total system is (may be in a mixed state due to topological ground state degeneracy) at zero temperature, and we compute ρ A , then there exists certain effective Hamiltonian on the EB such that The H EB is the effective Hamiltonian that causes the quantum fluctuations, and preserves the global symmetry broken from the gauge symmetries in A and B. The infinite temperature agrees with our conclusion that the region A is in equilibrium with maximal fluctuations on the EB. This may be related to the program of modular flows (see Ref. [38] and references therein).

Topological orders with gapped physical boundaries
So far we have studied topological orders on closed two-spaces. Nevertheless, in reality, materials fabricated to bear topological orders usually have PBs. In particular, non-chiral topological orders, which can be studied by lattice Hamiltonian models, have gapped PBs between the system and the vacuum. It is therefore important to extend our study to the cases with physical boundaries. To this end, we consider the extended QD model with gapped boundaries [29,30], which were developed to handle such cases. In this model, the bulk d.o.f. still take value in a finite gauge group G, while the d.o.f. on each PB take value in certain subgroup K of G. It is said that the PB is characterized by K. Our holonomy bases and fusion bases naturally extend to such cases, and the picture of quasiparticle fluctuations still applies. To be specific, in what follows, we consider two representative cases, seen in Fig. 7 and Fig. 8. The former is a cylinder with region A being a meridian strip, and the latter is a cylinder with region A being a longitudinal strip.

Cylinder case I
Consider the extended QD model on the cylinder in Fig. 7(a). The model is defined by a QD Hamiltonian in the bulk with gauge group G, and by a boundary Hamiltonian [30] on PBs specified by two subgroups K 1 and K 2 of G. The subgroup K 1 (K 2 ) dictates the gapped PB condition on the top (bottom) PB of the cylinder. Topological GSD occurs in general on a cylinder. As to be specified below, we will focus on only one particular ground state, which bears no non-contractible anyon loops. Similar to the case on a sphere, we choose a base point v A (v B ) in A (B). Since trivial holonomy is assumed everywhere, we only need to consider the d.o.f. on the EB, which is a disjointed union of two vertical line segments, each with two end vertices respectively on the top and bottom PBs. See Fig. 7(b) for an illustration. By a derivation similar to that for the ground state on a sphere, we can find one particular ground state (out of many others), which does not have any non-contractible loops, as where v labels the bulk vertices on the EB, and u = 1, 3 (u = 2, 4) label the vertices on the intersections of the upper (lower) PB and EBs. One can verify the K 1 , K 2 -gauge invariance at the PB vertices u and G-gauge invariance at v. Hence, |Φ cyl is indeed a ground state. Explicitly, the state is obtained by applying average (bulk and boundary) gauge transformations everywhere on the state with k u , g v set to be 1.
As in the sphere case, the ground state |Φ cyl is manifestly Schmidt decomposed; hence, the reduced density matrix is readily where K = K 1 ∩ K 2 . Similar as in the sphere case, one can show that ρ A is a projector multiplied by a normalization factor. Direct calculation leads to a normalized reduced density matrixρ where L is the total number of bulk vertices on the EB. The EE is See Appendix B for the detailed evaluation. Despite in this case the holonomy basis (4.1) has already made the EE computation simple, we would still like to write down the fusion basis as follows to manifest the picture of quasiparticle fluctuations.

The fusion basis
By examining Eq. (4.2) and (4.4), we can extract a picture of quasiparticle fluctuations. As illustrated in Fig. 7, The EB intersects PB at the four vertices labeled by u. The reduced density matrix consists of two projectors. One projector is which projects onto the subspace where we assign k 1 , k 3 ∈ K 1 on the upper PB and k 2 , k 4 ∈ K 2 on the lower PB. At other vertices on the EB, we assign group elements g v ∈ G. The other projector is an average of global right-multiplications: The reduced density matrix is a combination: This formula is similar to that in the sphere case, where the identity matrix in Eq. (3.13) is replaced by the projector P 0 . As suggested by the projector P 0 , a good fusion basis should be constructed in terms of the irreducible representations of K 1 , K 2 at the four vertices u and the irreducible representations of G at the other vertices v. The projector P G enforces a global right-multiplication symmetry, which leads to a fusion tree structure similar to Eq. (3.12). To be clear, we draw the fusion basis as where p 1 , p 3 (p 2 , p 4 ) label the irreducible representations of K 1 (K 2 ), and the j's are those of G. The x's and m's are the internal indices of the corresponding representations. Residing on the two big dots are the unitary transformations U j p,α , defined in Eq. (D.2), that decompose G representations to K 1,2 representations. On internal lines and trivalent vertices are the invariant tensor basis as defined in Eq. (A.3), for K 1 , K 2 , and G respectively.
It is more complicated to solve the eigen-problem of ρ A in such fusion basis than in the holonomy basis. One complexity is that such fusion basis is not orthogonal with respect to the d.o.f. at the two big dots (in terms of U j p,α ). (These dots are K 1,2 -morphisms but not G-morphisms, and hence are not invariant under P G .) We will not detail the computation of EE in such fusion basis.
The entanglement entropy in Eq. (4.5) consists of two terms, with the leading area term where L is the length of the EB in the bulk, i.e., the number of bulk vertices on the EB. The latter two terms reflects the contribution of quantum fluctuations at the four vertices at the intersections of the EBs and the PBs. Each such vertex on the top (bottom) PB allows |K 1 | (|K 2 |) types of quasiparticles. The subleading term is as a consequence of the global G-symmetry on all quasiparticles on the EB (including the four intersection points).

Cylinder case II
We consider another bipartite system on a cylinder as in Fig. 8(a). The extended QD model is the same as in the cylinder case I, where the bulk gauge group is G, and the top and bottom PB conditions are characterized by subgroups K 1 and K 2 of G. Again we will focus on only one particular ground state to be specified below. Configurations of group elements for a ground state. The dash lines are the two EBs. The group elements k 1 ∈ K 1 , k 2 ∈ K 2 , and g v , g ′ v , h ∈ G. We place the k's and g's in one region and their duplicates in the other, as a consequence of trivial holonomy condition.
Since region B now has two disjoint components, we will choose the base points in a way slightly different from that in the previous case. We pick two base points in region B, with v B1 on the top PB and v B2 on the bottom PB. See Fig. 8(a). The d.o.f. to specify a ground state are shown in Fig. 8(b).
The ground state of our concern takes the form The reduced density matrix is (4.14) which evaluates as (up to a normalization factor) In particular when G is Abelian, ρ A can be normalized to be a projector with trace being where K = K 1 ∩ K 2 . See Appendix C for the detailed computation. The EE in the Abelian case is thus which agrees with the result in Ref. [22]. In the non-Abelian case, we need some more tricks to obtain a closed-form formula of the EE. We leave the analytic computation in the fusion basis to the next subsection.

The fusion basis
We see that ρ A is a projector. We also see that ρ A is gauge invariant at all of the base points chosen, in kets and bras respectively. These facts imply all four conditions on the EB-states are fulfilled. In this subsection, we analyze the entanglement in the picture of quasiparticle fluctuations using the fusion basis. Let be the right-multiplication acting on g v and g ′ v by g 1 and g 2 . Define two projectors P G = 1 |G| g∈G R g,g (4.20) Then the reduced density matrix can be written as where P K 1 K 2 = P K 1 P K 2 . If we traced out region A, we would get The projectors impose the global constraints in the bulk and on the PB respectively. The P G projects onto the states with zero charge in bulk, while P K 1 and P K 2 projects onto the states with zero boundary k 1,2 charge at each boundary component.
We now try to find the eigenvectors of ρ A . First consider the effect of P G . Let us take a fusion basis, as obtained in Eq. (3.12): where we divide the labels j 1 m 1 , . . . , j L m L into two sets, with the d.o.f. on the upper (lower) EB all on the left (right) hand side of η 0 . Between these two parts there is one internal line labeled by η 0 , an irreducible representation of G. The combined projector P K 1 K 2 is a projection on η 0 and contribute a factor Which satisfies The eigenvalues are N η 0 , whose degeneracies are calculated in a way similar to that in Eq. (3.16), 27) whereη 0 in the sum means that η 0 is a fixed label and is not summed.
The trace evaluates as The normalized eigenvalues are then λ η 0 = N η 0 /W with degeneracy W η 0 , and hence the EE computes as (4.29) The EE again contains two terms, S E = S 0 + S 1 , with the leading area term S 0 = L log |G| and a subleading term as a consequence of three mutually un-commuting symmetries: P G acting on all quasiparticles and P K 1 (P K 2 ) acting on the quasiparticles on the upper (lower) EB. This closed-form formula is a novel result, which has not been obtained in other works on the EE of the QD model. In a recent work [22], a numeric result of the EE in this case for G = S 3 , K 1 = {1, 4}, and K 2 = {1} was given. Let us now verify as an example that our formula (4.1) yields the same numeric result. In this example, we have See Appendix E for S 3 and its representations. Eq. (4.1) quickly leads to which is precisely the result in Ref. [22]. When G is Abelian, N η is either 0 or 1 (as can be directly verified from the definition in Eq. (4.25)). Hence the last term in Eq. (4.30) vanishes, recovering the result (4.18).

Generic cases
We have discussed the picture of quasiparticle fluctuations in two representative cases on a cylinder; however, this picture is valid on more generic surfaces. We will summarize some typical features of the previous analysis and then extend them to a generic surface. For example, consider an open surface with PBs as illustrated in Fig. 4.3. We do not consider any genus in the bulk because we always choose a particular ground state such that the global (topological) d.o.f. will not affect our computation. Specifically, we choose a ground state in which the holonomy along all non-contractible loops are the unit element of G, which makes any genus invisible in the holonomy bases. Figure 9. (a) An open surface with boundaries. The thick lines are the EBs that separate the system into two regions A and B. The five PBs are specified by the subgroups K 1 , K 2 , K 3 , K 4 , and The subspaces with group elements k 1 , k 2 ∈ K 1 , k 3 , k 4 ∈ K 2 , and k 5 , k 6 ∈ K 3 , and G elements g v and g ′ v assigned to the EBs.
Let us try to construct the reduced density matrices on generic surfaces. We can see that the reduced density matrices (3.13), (4.8), and (4.22) all are a combination of three types of projectors. The first type includes projectors due to the boundary conditions on the PBs, such as the P 0 in Eq. (4.6), projecting all group elements at the intersection of PBs and EBs into the characterizing subgroups for PBs. (This type of projectors also depends on the choice of the ground state but we have already restricted to ground states without anyon loops throughout the paper.) The other two types are due to the gauge symmetries in regions A and B respectively, such as the P G and P K 1 K 2 in Eq. (4.22), where G is the gauge group in A, and K 1 (K 2 ) is the gauge group of the top (bottom) PBs in B.
These observations comply with the four conditions on the EB in Sec 3.3. Now we may formulate these conditions mathematically by writing down a generic form of the reduced density matrix.
Keep Fig. 9 in mind as an example, we itemize the rules to write down the reduced density matrix: 1. The first projector projects onto a subspace with only K 1 , K 2 and K 3 d.o.f. at all intersections between EBs and PBs, and assigns G-elements elsewhere. In this example, it is written as where u = 1, 2, . . . , 6.  (4.34) and

Suppose region
Since these projectors are due to the global symmetry in region B, we combine them as the projector 3. Repeat the second step also for all disconnected components of A. In the current example, we have (4.37) 4. We finally combine all the projectors and obtain This form of ρ A is a consequence of the four conditions on EB-states. Conditions 1, 2, and 3 allow equal-probable configurations of charge quasiparticles on the EB. Hence, we arrive at the P 0 in Eq. (4.3). Under condition 4, the two global symmetry constraints in A and B implies two projection operators P A in Eq. (4.37) and P B in Eq. (4.36). Combining all four conditions we arrive at the formula for ρ A in Eq. (4.38). Before explaining how we combine the projectors P 0 , P A and P B in such a way to get ρ A , let us figure out how the observers in A would view their world when B is traced out. They cannot distinguish a PB from an EB but think all the PBs and EBs are physical boundaries of their world. See Fig. 10. By Condition 4 on EB-states, we have a global symmetry acting on the EB d.o.f. due the broken gauge symmetry in B. The resulting global symmetry is implemented by P B . Observers in A cannot understand why they can observe such a global symmetry on a boundary that was an EB because they are blinded from B. In this way, the projector P B encodes how the topology and PBs of B may affect the EB-states observed in A.
There is another global symmetry acting on the EBs. Combining the already existing PB conditions (characterized by K 1 , K 2 , and K 3 ) and the effective boundary conditions (expressed as P B ) on the EBs, the gauge symmetry in A also breaks into a global symmetry, now imposed by P A . The above interpretation of P A and P B enables us to write down the combination of the projections in the order in Eq. (4.38).
We may choose to trace out region A instead and get in which the order of projectors accords with the above interpretation of the projectors.
Since P A and P B do not commute, ρ A = ρ B . The nonzero eigenspaces of ρ A and ρ B are not identical. This is reasonable because the regions A and B have different spatial topologies and different PB conditions. The local observers in A and B can distinguish such differences in the EB-states. For example, local observers will find the constraint (4.25) on the fusion basis (4.24) in Cylinder case II, and find the constraint on the dotted part in fusion basis (4.9). Nevertheless, the EE computed from ρ A and that from ρ B out to be equal.

The disk case
As a simple example, we apply our generic formula (4.38) to the disk case. We consider the extended QD model a disk as in Fig. 11(a). The bulk gauge group is G, and the PB is specified by a subgroup K of G. The system is bipartite by an EB (solid circle in the figure). The ground state on the disk is unique, as to be specified below. The d.o.f. would be g v for all v on the EB. The global symmetries would be P A = P G being the average of right-multiplications over G on g v for all v, and P B = P K over K. Hence the reduced density matrix ρ A is (4.40) In the second equality we have used P G P K = P K P G = P G . Hence the entanglement spectrum and entropy are the same as those in the sphere case.

Discussions
In this paper, we always choose a ground state in which there are no non-contractible anyon loops. In the sphere case, certainly there are no such loops anyway but on a cylinder, there may be because a topological order on a cylinder may have ground state degeneracy [30,39,40]. In general, we can choose an arbitrary ground state, with or without non-contractible anyon loops. In such cases, while the analysis and the main result will remain unaffected, we would need to slightly extends the four conditions on the entangled EB-states. Moreover, one may also consider topological orders with excited anyons, i.e., not in ground states. This would not affect the physics of EE, although special cares may be needed for defining the entangled EB-states. We shall leave such cases for future work.
In this paper, each EB accommodates only the quasiparticles that are pure charges. The reason is that each EB is composed of consecutive edges of the lattice; hence, the boundary condition is specified by G itself. This boundary condition corresponds to the smooth boundary condition on the PBs in the extended QD model [19,29,30]. Under the smooth boundary condition, the gapped boundary is a consequence of condensing all the flux quasiparticles on the PB, such that the gapped quasiparticle excitations on the PB are charges only. The correspondence between our EB and flux condensation on the PB of the extended QD model is more transparent to Alice, the local observer in region A, to whom the EB is effectively a PB. There are other choices of EBs, on which the quasiparticles are either pure fluxes, dyons, or of mixed types. According to Ref. [19], the TEE of topological orders is independent of the choice of EB. We thus do not consider other choices of EBs here. Moreover, in Ref. [32], the authors computed the EE between two different topological orders separating by a gapped domain wall by treating the domain as the EB. Since the gapped domain wall is a result of anyon condensation on the wall, the correspondence between our EBs and anyon condensation is natural.
Although our study was done for the (extended) QD model only, we believe our results can be extended to other models of topological phases, such as the twisted quantum double model [26,30], the Levin-Wen model [6,[40][41][42][43] in two dimensions, and the twisted gauge theory model in three dimensions [44,45]. In such extensions, however, the dimension of group representations should be extended to the quantum dimensions of the anyons. The latter is more general measure in the modular tensor category theories-the general algebraic theory of quasiparticles. We shall leave the generalization for future work.

A Fusion basis in the sphere case
We have used the fusion bases to represent the EB-states. Here we provide a concrete definition of the fusion bases. We define a fusion basis by

B Reduced density matrix in cylinder case I
We show that the ρ A (4.2) is a projector up to a normalization factor. We have where c runs over K = K 1 ∩ K 2 . To see the second equality above, we note that for example, which also constrains c ∈ K. We can evaluate the trace of the projection operator β −1 ρ A : where L is the total number of bulk vertices (those not on PB) of the EB. Hence the EE:

C Reduced density matrix in cylinder case II
We show that the ρ A (4.2) is a projector up to a normalization factor when the gauge group G is Abelian. We have ′ .
(C.1) When G is Abelian, we can define a projector with the constant factor β = We can evaluate the trace of the projector β −1 ρ A : where L is the total number of vertices on the EB. The EE in the Abelian case then is The group S 3 has three irreducible representations. A representative set of the irreducible representations is g = 1 g = 2 g = 4 ρ 1 (g) 1 1 1 The representation matrices for other group elements can be obtained by the group multiplications. In this setting, Eq. (4.31) was obtained by direct computation.