Entanglement Entropy of Topological Orders with Boundaries

In this paper we explore how non trivial boundary conditions could influence the entanglement entropy in a topological order in 2+1 dimensions. Specifically we consider the special class of topological orders describable by the quantum double. We will find very interesting dependence of the entanglement entropy on the boundary conditions particularly when the system is non-Abelian. Along the way, we demonstrate a streamlined procedure to compute the entanglement entropy, which is particularly efficient when dealing with systems with boundaries. We also show how this method efficiently reproduces all the known results in the presence of anyonic excitations.


Introduction
Entanglement entropy has been a powerful probe to detect properties of matter. For example, entanglement entropy can probe different topological orders. In 2+1 dimensions, the entanglement entropy of a topologically ordered quantum system, when the system is divided into two regions A and B, takes the following generic form where L A is the length of the boundary of region A, the UV cutoff, and γ, the universal term named "topological entanglement entropy" [1,2]. It is known that when S(A) is evaluated on the ground state of the system on a sphere, for some region A consisting of n disconnected disks, and D the total quantum dimension of the topological order, defined as and d i are the quantum dimensions of the anyons a i of the topological order.
In the current paper, we shall inspect the entanglement entropy in a topologically ordered system (to be called a topological order for simplicity unless otherwise stated) with non-trivial boundary conditions.
To avoid any confusion, we shall refer to a boundary between two regions in a system as an entanglement boundary (EB) and a physical, gapped boundary between the system and vacuum a physical boundary (PB).
Non-chiral topological orders can admit gapped, or alternatively topological boundary conditions. In 2+1 dimensions, each such boundary condition is related to an algebra often called the "Lagrangian subalgebra" [3,4,5,6]. It is also characterized by a Frobenius algebra [7] in the modular tensor category describing the topological order or by a Frobenius algebra [8,9,10] in the unitary fusion category describing the fundamental degrees of freedom of the topological order. These boundary conditions are also known to correspond to the physics of "anyon condensation" [11,12,13,14]. Each such boundary is also associated to some modular invariant [3,7,15]. We would like to inspect whether the entanglement entropy can detect these boundary conditions, and if so, whether the resultant values correspond to certain topological invariants.
We will inspect this problem in the context of the twisted quantum double (TQD) models of (2 + 1)d topological orders [16], which are a Hamiltonian extension of the Dijkgraaf-Witten topological gauge theories [17]. The entanglement entropy of these models have been computed before, notably in [18,19] in the absence of physical boundaries. In the case of Z 2 toric code model, the case with PB was also briefly discussed in [20]. We will revisit the problem, and introduce a streamlined method. It involves reducing the problem systematically to one that is independent of system size. We also clarify the construction of Schmidt decomposition by systematically choosing a convenient canonical set of basis. As we will see, the modified discussion would enable an efficient and clear inspection of the scenario when we have non-trivial boundary conditions.

TQD (Dijkgraaf-Witten) models
Dijkgraaf-Witten topological gauge theories were formulated initially as a way of defining a path-integral of a discrete version of Chern-Simons theory. They were later adopted for defining quantum Hamiltonians, whose ground states admit exotic properties, that we now identify as the fixed-point wavefunction of topological orders. For a detailed discussion, we refer the readers to [21,16].
We only collect the necessary ingredients of the TQD models that would facilitate our exposition in the following. We note that the model is defined on a lattice Γ. The lattice does not have to be regular but without loss of generality, we shall consider a square lattice. There is a Hilbert space H l defined at each link l.
We first consider a special subset of the TQD models, where there is no twist. Such models are called the Kitaev models or quantum double (QD) models. A QD model with a finite gauge group G, dim H l = |G|, where one can choose a basis such that each basis state can be labeled by a group element of the group G, has the Hamiltonian (2.1) The subscript v denotes vertices, and p denotes plaquettes. The action of A v (g) and B p (g) is illustrated in figure 1.
Since all the operators A v and B p commute, the ground state can be generated by where |Ω is some appropriate reference state satisfying B p |Ω = |Ω . When the (closed) two dimensional space on which the state lives is a sphere, the reference state |Ω is simply given by all links taking the state corresponding to the identity element e of the group G, Figure 1: An illustration of the action of A v (g) and B p (g).
where l runs over all links of the lattice Γ. This state |Ω is of course a direct product state. When the 2-dimensional space has non contractible cycles, there would be a ground state degeneracy, and a basis of the degenerate ground states can be constructed by taking a set of reference states with closed ribbon operators acting on non-contractible cycle in the trivial reference state |Ω .
The analysis in the following for individual such states are all the same. We will discuss general linear combinations of these basis in later sections when we encounter the geometry of a cylinder.

Gapped boundaries
The PB conditions of the QD models have been discussed in [22,23] and subsequently generalized to the TQD models [9,24] , to the Levin-Wen models [8,10], and to higher dimensions [25].
We will illustrate our methods mainly using the QD model. Each PB is characterized by a subgroup K ⊂ G [23]. In fact, it is shown in [9] that even for a TQD model with a gauge group G, a PB condition is fully characterized by a subgroup K ⊆ G.

The PB Hamiltonian is given by
where A v B acts on the links connected to the vertex located at a PB given by , and B l B is a projector on the PB links to the subgroup K.
A ground state is generated in a similar manner as in (2.2), except that for vertices on the PB we replace a generic A v by A v B .

Revisiting the entanglement entropy of the QD models
In this section, we revisit the problem of computing the entanglement entropy in the QD models. We shall lay out a procedure that is improved compared with the discussion in [18,19,26]. The procedure would allow one to obtain the entanglement entropy in the presence of PBs in a systematic and clear manner.
Now for simplicity, we consider again the case of an entangling region R taking the shape of a disk on the sphere.
On the sphere, any topological order has unit ground state degeneracy. The ground state on the sphere is generated as described above, in equation (2.2). It is known that the operators A v acting on v away from the EB between the region R and its complement R do not contribute to generating entanglement between the regions. It is also known that the entanglement arises from operators A v b that act on the vertices v b along the EB and hence affect both region R andR at the same time. Hence, one can simply label some i-th EB configuration by the set of vertex operators Consequently, one can just focus on different EB configurations {A v b (g i b )} in each term in the ground-state wavefunction. As usual, there is a physical ambiguity over the definition of entanglement entropy on a lattice gauge theory. But here, we have taken the viewpoint explained for example in [27,28,29,30,31], and work with the extended Hilbert space, which should agree with the electric center in terms of the choice of operator algebra in the original gauge theory [32].
The game is to obtain a Schimidt decomposition to recover the reduced density matrix. Given that only A v b are responsible for the entanglement between R andR, a naive Schimidt decomposition is obtained as where L is the length, i.e., the number of links, along the EB. We note that the rhs above is indeed a direct product state in R andR, agreeing with the expression on the lhs, since b=1 is a tensor product of operators acting on R andR respectively. Had each EB configuration b A v b (g b ) always led to a pair of orthogonal internal states |R i and |R i , the Schidmit decomposition would have been completed. This is however not the case. Indeed, two different sets generating states |R i,j and similarly |R i,j for different pairs of i, j are generically not linearly independent. That is, our naive labelling of the internal states based on the EB configuration of A v b does not in general lead to |G| L orthogonal |R i and |R i separately.
There is a complication, namely that many of the |R i (|R i ) for different i correspond to the same state. This is because the global action of the internal A v∈R (R) in each connected As we are going to see, more generically when the region R orR has topologies more non-trivial than a disk, or when the system is placed on surfaces beyond a sphere, the ground state is still generated by in an analogous manner as in 2.2 with the reference state |Ω potentially dressed by ribbon operators (the discussion of these dressed states are postponed until section 4.6). The paremetrization adopted in 3.1 and 3.2 can still work quite generally. Generically, we would be met with a new complication. That is, whenever we transform a state by A v (g) at all v within region R, including the entanglement boundary v b at the same time, we would keep |R i invariant, while taking |R i to another state |R j . The EB configuration after such a transformation has been shifted This implies that some of the |R i,j with i = j may in fact be the same but they do connect to different |R i,j , leading to entanglement pattern in the wavefunction of the form |R i ⊗ (|R i + |R j + · · · ).
In the classic literature on the subject [18], further analysis is based on defining a huge group G corresponding to the action of A v (g) on the entire lattice, and attempting to obtain the quotient group G R,R , where we quotient by the action of the corresponding groups G R and GR, where G X includes action of A v on vertices v within X, excluding the EB. This is of course Mathematically correct but the analysis on complicated situations where there are PBs can be confusing at times.

A modified analysis
The improvement proposed in the current paper is to consider explicitly a division of the collection of |G| L EB configurations b A v b (g b ) into distinct sets. The choice of such a division is such that each group would contribute to identical blocks in the reduced density matrix. i.e. The complications that arise due to global applications of A v (g) within R or R described above could generate off-diagonal elements only within each block.
These independent blocks are obtained as follows. Let us begin with the simplest scenario, where the model is defined on a 2-sphere. Region R is a disk with circumference L on the sphere. That is, the EB between R andR consists of L links. This is illustrated in figure 2.
Then there are |G| L−1 sets. Each set is obtained in this way: take a configuration along the EB as a reference, which is a result of acting the vertex operators Π i∈EB A v i (g i ), then multiply each g i with one and the same group element g ∈ G to reach another configuration in the set, and repeat this for all element of G. This is the G 1 -orbit of configurations along the EB and contains precisely |G| configurations referred to above. Each configuration in a G-orbit is clearly labeled by certain group element g ∈ G. According to the discussion above therefore, we write where g i labels the global shift along v b over a reference representative EB configuration in a G-orbit. Now in this case, the entanglement boundary is contractible both within R, and inR. We note that each extra global action of g i at A v b is to create a pair of g i shifts that form a closed loop in both R andR simultaneously, and they are contractible since the EB is contractible. This immediately suggests that |R (R) (g i ) = |R (R) (g j ) for all i, j, due to the separate action of vr∈R (R) A vr within the regions.
We thus conclude that the wavefunction within this G-orbit can be simplified to |ψ 1 block = |R(g 1 ) ⊗ |R(g 1 ), (3.4) which is a direct product state within the block. Therefore each block contributes to a 1 dimensional projector to the reduced density matrix.
The entanglement entropy thus reads which is just the log of the number of blocks and recovers the well-known result.
Now, consider an EB of N disconnected components, each component contractible in both R andR. For any given configuration of i.e. we collect EB configurations related to each other by at most a global shift by some element g i on each individual EB. There are thus |G| N members in a set, which is a disjoint union of N G-orbits. There are |G| L−N separate G N -orbits that would not interfere with each other as we mod out actions of global transformations in various regions. The reduced density matrix would be reduced to a block diagonal form with |G| L−N blocks, each being identical.
Individual member within a G N -orbit can thus be labelled by an N -tuple: (g 1 , · · · g N ).
Using the same reasoning, for N disconnected EBs, we have |G| L−N blocks, while each block is again equivalent to a single direct product state. We thus obtain which is indeed the well known result for the entanglement entropy.
This method can also be readily applied to compute entanglement entropy in the presence of anyon excitations. We will for illustrative purpose demonstrate these applications in the appendix.
To summarize, this method systematically reduces a problem of treating a huge density matrix that scales as |G| L for an entanglement boundary of size L, to one whose dimension only scales at most with |G| N , where N is the number of disconnected components of the entanglement boundary, allowing a clear analysis even in complicated situations.
We will now apply this set of methods to the case where there are PBs and where the analysis can get substantially more complicated and the advantage of the method more pronounced.

Entanglement entropy for different PBs
In this section, we would like to apply our trick of entanglement entropy computation to the case where there are PBs. The choice of the PB condition would modify the final results of the entanglement entropy. As to be seen, the topological entanglement entropy has a subtle dependence on the PBs.

Case I: Region R is a disk away from the PB of a disk
Consider a disk whose PB is characterized by a subgroup K ⊆ G. Then consider a region R that is also a disk away from the PB. This is illustrated in figure 3 for R consisting of a single connected component.

Region R
Region ഥ R Figure 3: a disk away from the PB of a disk In this case, we have one EB, so we divide the EB configurations into |G| L−1 G-orbits, each G-orbit has |G| members labeled by (g). Now, for precisely the same reason as in the previous analysis of a region with a disk topology, all |R(g i ) = |R(g j ) for all i, j. This immediately implies that within each G-orbit we have a direct product state.
The entanglement entropy is thus still given by counting the number of individual blocks, and recovers the result We have generalized the result directly to the case where R contains N disconnected disks. The result is insensitive to the presence of the non-trivial gapped boundaries.
We note however, that a new ingredient has crept in in the current situation. Although it did not change the entanglement spectrum, it would make a difference in later analysis. The new ingredient is that the EB may not be contractible inR. Therefore, not all |R(g i ) are the same. It is interesting to check which of the |R(g i ) are in fact orthogonal. We note that a global action of A v∈R (k), k ∈ K can generate a closed loop of k shifts along the links connecting to the EB. This means that for all k ∈ K. i.e. All |R(g i ) for g i belonging to the same left coset of K corresponds to the same state. For completeness, we thus have where c i is a representative of a left coset. There are |G|/|K| left cosets.

Case II: Region R touching the PB of a disk
Here we continue to keep the state on a disk with a PB characterized by K ∈ G. The region R however, touches the PB. The EB is thus a line that begins and ends at the physical boundary. This is illustrated in figure 4.

Region R
Region ഥ R Figure 4: Region R touching the PB of a disk In this case, if the EB contains L vertices, then two of them sits at the PB, while the rest are located in the bulk. Therefore, the total number of configurations of possible A v sitting at the EB is given by |G| L−2 |K| 2 . Naively we would like to divide these configurations into G-orbits. Nevertheless, we are not able to do so here because A v B (g) located at the PB are restricted to g ∈ K. Instead, we can divide these configurations into K orbits, with |K| members in the orbit, each labeled by (k). There are thus |G| L−2 |K| distinct orbits.
In this case, the EB is not contractible in either R orR.
From the analysis of the previous section, we note however that not all the internal |R(g i ) or |R(g i ) are independent. They again satisfy for all k ∈ K, and thus all |R One therefore concludes that the orbit is contributing to one direct product state.
The entanglement entropy is then given by where the first term is grouped together and taken as the area term, and the second term a topological term resulting from a change in the global constraints. We see the first indication that the topological entanglement entropy is indeed sensitive to the physical boundary conditions.

Case III: Region R being a vertical slit on a cylinder with two PBs
Now consider a state on a cylinder with two PBs characterized by two subgroups K 1 and K 2 of G respectively (see Fig. 5).
Now on a cylinder with a non-contractible cycle, generically there would be degenerate ground states and the entanglement entropy would depend on the precise linear combination of ground states.
It is known that one can construct a set of basis states in the degenerate groundstate subspace using ribbon operators. This is obtained by wrapping ribbon operators around non-contractible cycles. In the case where there are PBs, one can also construct basis by attaching ribbon operators that stretch between the upper and lower PBs. Not all ribbon operators can end at the PB without leading to boundary excitations, except those corresponding to condensed anyons at the PB.
We first take the simplest basis state, corresponding to no ribbon, in which the state is still describable by equation (2.2). We will postpone the discussion of generic ground states to section 4.6.
Region R Region ഥ R Figure 5: A vertical slit on a cylinder The first scenario we consider is a region R that corresponds to a strip that connects the upper PB K 1 and the lower PB K 2 . In this case, there are two disconnected EBs too. There are altogether |G| L−4 |K 1 | 2 |K 2 | 2 different EB configurations. For similar reasons as in the previous subsection, we are not able to divide these EBs into G 2 orbits because of the restriction of the physical boundary vertex. We are however allowed to divide the EB into K 2 -orbit, where K is the intersection of K 1 and K 2 , which is itself a subgroup of G. Each member in the orbit is now labeled by (k 1 , k 2 ). In this case, each connected component of the EB is not contractible, either in R orR.
As in the previous examples, we now systematically proceed in two steps. First we deterimine which of the |R (R) (k 1 , k 2 ) are identified. Then we inspect global actions in each connected component of R orR to look for potential off-diagonal terms in the Schimdt decomposition.
From similar analysis of non-contractible EB in the previous subsection, we conclude that where K are group elements shared by K 1 and K 2 .
Then we look for global actions allowed within R orR that keep the respective region invariant. Again that is restricted to elements k in K, which also happens to take |R (R) (k 1 , k 2 ) → |R (R) (k 1 k, k 2 k) . As a result, following our procedure in the previous examples, we conclude that each K 2 -orbit breaks up into |K| entangled states. The entanglement entropy takes the form For region R made up of N such strips connecting the top PB to the bottom PB, there are 2N disconnected EBs. Repeating exactly the same analysis, we find that the entanglement entropy would take the general form Now we consider another interesting case. The state is still defined on a cylinder whose PBs are characterized by groups K 1 and K 2 respectively.
However the region R is now taken as a strip wrapping the non-contractible cycle of the cylinder, separated from both PBs.
Again, we will begin with an analysis based on the state (2.2). It turns out that this case is one that is the most interesting we have encountered so far.
The EB is made up of two disconnected components each non-contractible in either R andR. As before, we first divide all the EB configurations into G 2 -orbits with each member labelled by (g 1 , g 2 ), where g 1 corresponds to the upper EB and g 2 corresponds to the lower EB. The total number of EB configurations is |G| L , where L is the number of vertices on the EB, so there are |G| L−2 G 2 -orbits and each G 2 -orbit has |G| 2 members.
Again, our first step is to determine which of the |R (R) (g 1 , g 2 ) are dependent. From action within R one can generate a pair of g ribbons with A v (g) near the upper and lower entanglement boundary simultaneously. Similarly A v in regionR can generate a ribbon k ∈ K 1 in the upper entanglement boundary and k ∈ K 2 in the lower entanglement boundary. Therefore we have |R(g 1 , g 2 ) = |R(g 1 g, g 2 g) , |R(g 1 , g 2 ) = |R(g 1 k 1 , g 2 k 2 ) . (4.10) Then we would like to analyze how these independent basis states of R andR are entangled. As in the previous analysis, global actions in R together with the EB vertices preserve R but shiftR. It takes the EB configuration from (g 1 , g 2 ) → (g 1 g, g 2 g), thus taking the states |R(g 1 , g 2 ) → |R(g 1 g, g 2 g) . This means that |R(g 1 , g 2 ) is paired with both |R(g 1 , g 2 ) and also |R(g 1 g, g 2 g) . These global actions form a group G. Similarly, a global action inR and EB preservesR while taking the EB configuration (g 1 , g 2 ) → (g 1 k 1 , g 2 k 2 ), and thus taking |R(g 1 , g 2 ) → |R(g 1 k 1 , g 2 k 2 ) . These global actions form a group K 1 ⊗ K 2 .
To recover a Schimdt decomposition and subsequently the reduced density matrix, we need to count the number of times each independent |R(g 1 , g 2 ) gets paired with an independent |R(g 1 , g 2 ) after taking into account the redundancy (4.10). To that end, we further divide these |G| 2 EB configurations (g 1 , g 2 ) into sets of |G|. i.e. The members in each G 2 -orbit is further divided into |G| different G-orbits. Members in the G-orbit can be related by (g 1 , g 2 ) = (g 1 g, g 2 g). The EB configurations (g 1 , g 2 ) and (1, g 2 g −1 1 ) are in the same G-orbit. Therefore we can take (1, r) where r runs through the group G as the representative of a G-orbit in each G 2 -orbit. Other members of the G-orbit containing the representative (1, r) can be parametrized by (g, rg).
From (4.10), it implies that all members in a given G orbit are attached to the same |R(g 1 , g 2 ) state.
Next, we would like to analyze the effect of those global actions K 1 ⊗ K 2 that preservē R on the EB configurations. Under these actions, members in each G-orbit are allocated to the other G-orbits. For simplicity, consider the G-orbit represented by (1, 1). The members in that G-orbit are denoted by (g, g). For a specific choice of k 1 and k 2 , (g, g) is mapped to (gk 1 , gk 2 ). This is in the G-orbit represented by (1, gk 2 k −1 1 g −1 ). ie. There exists ag such that (gk 1g , gk 2g ) = (1, gk 2 k −1 1 g −1 ). (4.11) That is to say, under the action of k 1 and k 2 , members in the G-orbit (1, 1) are mapped into G-orbits labelled by the elements in the conjugacy class of k 2 k −1 1 . The number of members mapped into each G-orbit is equal to the order of the centralizer of k 2 k −1 1 : |Z(k 2 k −1 1 )|. Since K 1 ⊗ K 2 are actions that preserveR, it implies that |R(1, 1) and |R(1, gk 2 k −1 1 g −1 ) are paired with the same state |R (1, 1) . The analysis can be carried out for other G-orbits by the replacement (1, 1) → (1, r). For another pair (k 1 , k 2 ) in the group K 1 ⊗ K 2 satisfying k 2 k −1 1 = k 2 k −1 1 , the re-distribution of members of a G-orbit would be identical to that resulting from the action of (k 1 , k 2 ). If k 2 k −1 1 and k 2 k −1 1 belong to different conjugacy classes, (k 1 , k 2 ) will map the members in the G-orbit (1, 1) into G-orbits labelled by elements in a different conjugacy class other than the conjugacy class of k 2 k −1 1 .
Consider also the case where k 2 k −1 −1 for someg ∈ G. From the analysis above, it implies that members of a G-orbit would get mapped to the same destination G-orbits under the actions of (k 1 , k 2 ) and (k 1 , k 2 ). What is of note is that while these actions share the same destination G-orbit, the actual collection of destination members for the two actions have no overlap. That is, if the set K 2 K 1 -generated by all pairs k 2 k −1 1 for k i ∈ K i -contains N elements belonging to the conjugacy class C, the number of members in the G-orbit (1, 1) mapped into the G-orbit labelled by an element of that conjugacy class would be N · |Z C |, where Z C denote the centralizer of a representative element in the conjugacy class C.
In the following, we would like to prove this assertion. Suppose {c i } are the elements of a conjugacy class C of the group G. Denote the centralizer of each c i by Z(c i ). For every pair of c i and c j , choose a specific q i,j ∈ G such that c i = q i,j c j q −1 i,j and q i,j = q −1 j,i . Specifically, let q i,i = 1. Then for each c i , the elements in the set q j,i Z(c i ) would map it to c j . Therefore, we only need to show that for all the different i, the sets q j,i Z(c i ) have no overlap.
We observe that Z(c i ) = q i,j Z(c j )q −1 i,j . So the set q j,i Z(c i ) can be written as Z(c j )q −1 i,j . For different i the sets Z(c j )q −1 i,j are the right cosets of Z(c j ), so the sets Z(c j )q −1 i,j form a partition of the group G, completing the proof.
We note that there is a redundancy in k when K 2 ∩ K 1 contains more than the identity element. Collecting the observations above, we are ready to write down the reduced density matrix of region R. Suppose the set K 2 K 1 contains N elements belonging to the conjugacy class C, then N · |Z C | members of the G-orbit (1, 1) are mapped into each G-orbit labelled by the elements in the conjugacy class C.
The block of the reduced density matrix is (ρ R ) 1 block = trR g,g ,r,r ∈G |R(g, rg) |R(g, rg) R(g , r g )| R (g , r g )| = g,g ,r,r ∈G R (g , r g )|R(g, rg) |R(1, r) R(1, r )|. (4.12) Thus the coefficient of the term |R(1, r) R(1, r )| is g,g ∈G R (g , r g )|R(g, rg) . From the analysis above, the summation g,g ∈G R (g , r g )|R(g, rg) should equal to N · |Z C | multiplied by the dimension of the set K 1 ∩ K 2 . Since every entry of the reduced density matrix has a factor of the dimension of the set K 1 ∩K 2 , it will disappear after normalization and doesn't affect the final value of the entanglement entropy. Therefore, the procedure of writing down the reduced density matrix can be summarised as following. We only need to count the number of elements of K 2 K 1 belonging to each conjugacy class of G. Then for each conjugacy class C we just put the number N · |Z C | on the proper places in the first line in the block of the reduced density matrix. The analysis for other G-orbits is similar and thus we can write down the block of the reduced density matrix line by line. Other lines are just some permutations of the first line. Recall that the full reduced density matirx is obtained by |G| L−2 such blocks, we can calculate all of its eigenvalues within the G 2 -orbit.

Abelian groups
For Abelian groups, the g appearing in the r.h.s of (4.11) clearly cancels. Hence, the entire G-orbit (g, g) is mapped to one and the same G-orbit labeled by (1, k 2 k −1 1 ). This implies that the wavefunction takes the following form: The group element k however runs over the group generated by group multiplications of the elements of K 1 and K 2 . If K 1 ⊆ K 2 then this gives K 2 . Otherwise it is denoted more generally by K 1 K 2 already introduced in the discussion above. (We note that while in a non-Abelian group G, K 1 K 2 does not generically form a group, it does in the current case of an Abelian group G. ). The product kr i as k varies over K sweeps through the right coset of K in G. Thus the sum over r i only runs from i = 1 → |G|/|K 1 K 2 |, picking a representative of each distinct coset of |K 1 K 2 |.
We thus conclude that for the case of Abelian groups, the entanglement entropy of this region reads (4.14)

Case V: Multiple PBs
Region R K 1 K 2 K 3 If there are multiple PBs, the computation of entanglement entropy will be more complicated. Here we introduce how to compute the entanglement entropy in the case of three PBs. Suppose the PBs are characterized by subgroups K 1 , K 2 , K 3 ⊆ G, and suppose the region R is chosen as in Fig. 7. Thus the EB has three disconnected components We divide the EB configurations into G 3 -orbits with each member labelled by (g 1 , g 2 , g 3 ). Global shifts inR take (g 1 , g 2 , g 3 ) into (g 1 k 1 , g 2 k 2 , g 3 k 3 ) and global shifts in R take (g 1 , g 2 , g 3 ) into (g 1 g, g 2 g, g 3 g).
First we divide a G 3 -orbit into |G| 2 subsets according to the global shifts in R. Then each member in a subset can be labelled by (g, r i g, s j g), where r i and s j run through the group G. The action of k 1 , k 2 , k 3 takes (g, r i g, s j g) into (gk 1 , r i gk 2 , s j gk 3 ). Right multiplying the three terms by k −1 1 g −1 , we get (1, r i gk 2 k −1 1 g −1 , s j gk 3 k −1 1 g −1 ). As g runs through the group G, gk 2 k −1 1 g −1 and gk 3 k −1 1 g −1 run through the conjugacy class of k 2 k −1 1 and k 3 k −1 1 . Due to the lack of a general relation between the centralizers of k 2 k −1 1 and k 3 k −1 1 , we can't deduce a general procedure to write down the reduced density matrix. But for a specific group G, one can still write down the reduced density matrix and compute the entanglement entropy.

More generic basis states with non-trivial wrapping ribbons
In the above analysis, we have considered only the entanglement entropy of a specific reference ground state (2.2, 2.3). On a manifold with non-contractible cycles, or open boundaries where ribbons could end, the ground state becomes degenerate. In the case of a cylinder, we can construct a complete basis of the ground states by wrapping magnetic ribbons around the non-contractible cycles, in addition to magnetic ribbons in the axial direction connecting the top and bottom physical boundary on the trivial state |Ω in (2.3). The ribbon operators are discussed in the appendix figure 11. In the case of pure magnetic ribbons, we simply sum overall g in the projector with equal weights. The end result is simply that we give up the projection. Since we are acting the ribbon on the reference state (2.3), it amounts to a shift of a set of links by some group element h. Such a ribbon can wrap the non-contractible cycle, in which case h ∈ G. When the ribbon is one that connects the upper and lower boundaries of the cylinder, h ∈ K 1 ∩ K 2 . We can of course have independent ribbons wrapping the non-contractible cycle and extending in the axial direction. We can thus label a reference state by a pair of group elements (h e , k a ), where h e denotes a ribbon wrapping around the non-contractible cycle, while k a extends along the axial direction. In the presence of both, we have further restrictions of h e and k a . Namely, h e and k a should commute.
A generic ground state basis state is thus obtained by where |Ω(h e , k a ) corresponds to a reference state obtained by action of the corresponding ribbons on the trivial reference state (2.3). There is one extra complication in the presence of physical boundaries. These states are not all independent. In the presence of physical boundaries, as we have seen, a global action of v B A v B (h) essentially generates a noncontractible ribbon corresponding to the group element h. Therefore, in the presence of these projectors in (4.17), we have for k i ∈ K i , where K 1 is the subgroup characterizing the top physical boundary, and K 2 the bottom physical boundary. These redundancy has already been discussed in [14] which analyzes the number of degenerate ground states in the presence of physical boundaries.
To compute entanglement entropy, we note immediately that the result for each individual reference state is completely independent of the particular choice of these basis state.
A generic linear combination of these states, however, merits extra analysis. Since the situation in section 4.4 is the most interesting one, we will take it as an illustration.
To simplify the discussion further, we will restrict our explicit examples to Abelian groups. For given G and K 1,2 , we first construct the ground states. The set generated by K 1 and K 2 was denoted K 1 K 2 in the previous section. Here, we will call it K 1 K 2 = K which is also a subgroup of Abelian G. The intersection is denoted X ≡ K 1 ∩ K 2 . The ground state subspace is thus |G|/|K| · |X| dimensions. From the discussion above, a generic state is given by where r are group elements of the quotient group G/K, which can be treated as a representative of the coset of K. The second label x denotes axial ribbon and x ∈ X.
Ribbon x cuts through all the regions involved, while ribbon r can lie in any of the regions. The initial position of r does not matter, since the product of projectors A v serves to deform it in all possible topologically trivial way. Therefore, without loss of generality, we can directly choose to pick the reference state with r residing within region R.
For each individual basis, we first analyze it exactly as we did in the previous sections, dividing each |Ψ(r, x) into linear combinations of G 2 orbits, and organize our states into the form as in (4.13) where we keep track of (r, x) in the subscript of the states. The most important ingredient in the remaining analysis is that these states |R (R) (g 1 , g 2 ) (r,x) may not be independent for different (r, x). There is a set of simple relations between them. There is a correspondence of states between states in a given G 2 orbit. Since we have used the ansatz where the ribbon r resides in region R in the reference state, it implies the following relation |R(1, k p r i ) (r 1 r 2 ,x) = |R(1, k p r i r 1 ) (r 2 ,x) , |R(g, gr i ) (r 1 ,x) = |R(g, gr i ) (r 2 ,x) , (4.21) for all r i ∈ K. Note that states in different x sectors are immediately different and orthogonal, following from topological reasons -that the axial ribbon x always penetrates each region an even number of times i.e. whatever goes in comes out.
Then we would like to obtain the reduced density matrix. Tracing outR, it gives As we already anticipated above, the reduced density matrix is diagonal in x.
It only remains to analyze each term for x fixed. Using (4.20), we then have where we have used (4.21), and simplified notation by replacing Here, r i used to denote a representative in the coset of K in (4.20) and it is now simply denoting a group element of the quotient group G/K. One can immediately see that the state with maximal entanglement would be the one where c r 1 = 1 while all other c r j = 0. Each of the |G|/|K| state would contribute to entanglement, and we recover our previous result (4.14). The minimally entangled state would correspond to having all c r i equal. Then within the G 2 orbit we have a direct product state. In which case the entanglement entropy is given by As discussed in [20] this is a generic method of recovering the "anyon line" eigenbasis. We note that there is now non-trivial dependence on the boundary conditions for the maximally entangled state. The result gets more complicated without a clean closed form expression in the case of non-Abelian theories, although the minimally entangled shares exactly the same entanglement as in (4.25).
A complete understanding of the physical interpretation of these numbers would require new examples involving more general topological orders beyond lattice gauge theories.

A physical interpretation
In the above, we explained how to recover the Schmidt decomposition by dividing EB configurations of A v into distinct orbits or sets, and with members of each set labeled by group elements g. We would like to comment on the meaning of these labelling. Within the same set, different configurations are related to each other by an overall action of A v (g) for any g ∈ G along each given connected entanglement boundary component. The action of A v (g) on a closed loop is equivalent to creating a pair of closed (magnetic) ribbons of type g, one inside region R and the other in regionR. Therefore each member in the set really is a set of entangled state, for a single connected EB where |Ω is some reference state, the subscript i denotes the particular connected component among all components of the EB, and |Ω is some reference state (in this case it is the ground state). These states |g i are to denote the creation of extra magnetic ribbon in region R on top of the reference state. Now, had these |g i been orthogonal to each other, the above expression is a Schmidt decomposition. The entire analysis discussed in the current paper is about deciding which of these |R(R) (g i ) are in fact orthogonal to each other given that the reference state is a ground state generated by all possible linear combinations of A v covered on top of the direct product state. In the case of non-trivial topologies we have to include ground state basis states built from (magnetic) ribbons wrapping non-contractible cycles which are then buried under all possible A v showered on top. But the bottom line is that the summation of A v makes some of these |g i R(R) linearly dependent.
The key piece of physics is that for all other A v acting on the interior of R orR, they necessarily create ribbons that are contractible and hence topologically trivial. Therefore creation of g i in R(R) can generically be undone by A v from within the respective region if the g i ribbon (and thus the particular component of the entanglement boundary) is topologically trivial. This is the gist of the physics of the topological entanglement entropy. It is determined by the Gauss constraint which can be reformulated as having a bath of magnetic ribbons rendering these g i states linearly dependent.
What is interesting in the analysis above, where there are PBs, is that the content of this bath of magnetic ribbons get modified. By restricting to a subgroup K sitting at the PB, it is physically equivalent to allowing some ribbons to be created and destroyed at the boundary. This is of course a manifestation of anyon condensation at the boundary [22,14]. This changes the bath of ribbons in R(R) depending on their orientation in relation to the physical boundary, ultimately modifying the topological entanglement entropy. The interplay of leaking of ribbons at the boundary and the bulk bath of ribbon is clearly visible in the analysis. It appears that the interplay displays an interesting and complicated pattern when it comes to a generic non-Abelian group with multiple disconnected PBs and EBs. This issue has been pointed out before in [14], where they built ground state basis in the presence of multiple boundaries characterized by different anyons condensation connected by a bulk. Then, the treatment was by no means systematic, and it was dealt with in a case-by-case manner, based on the precise mutual statistics of the condensed anyons and their fusion properties. The results in the current paper may shed new light to this problem, leading to a complete understanding. Next consider the case that the lattice has a physical boundary with non-trivial cocycles ϕ : K × K → U (1). It satisfies the 2-cocycle condition:ϕ(kl, m)ϕ(k, l) = ϕ(k, lm)ϕ(l, m) [23]. The result is the same and it can be shown that the cocycles doesn't affect the computation. First consider two vertexes 1, 2. Suppose the A v operators on these two vertexes are A 1 (a), A 2 (b), originally. Then the values on the edges are a, ab −1 , b −1 , c correspondingly, and we have a phase factor ϕ(a, b −1 ) −1 . Next we replace the A 1 (a), A 2 (b) operators by A 1 (ag), A 2 (bg) as above, where g is an element of the subgroup K. We apply A 2 (bg) first and then A 1 (ag). This is illustrated in figure 9. The values on the edges become ag, ab −1 , g −1 b −1 , and the phase factor becomes ϕ(ag, g −1 b −1 ) −1 . We have ϕ(a, g)ϕ(ag, g −1 b −1 ) = ϕ(a, b −1 )ϕ(g, g −1 b −1 ). (4.27) We also have ϕ(g,

Comments on the twisted boundaries
The new phase factor differs from the original one by a factor ϕ(b, g)/ϕ(a, g). If we do the above for all the vertices on the physical boundary of region A, most terms will cancel and only two terms which are independent of the vertices inside region A are left. So we can still use the same method to compute the entanglement entropy as if there's no non-trivial cocycle on the boundary. It is observed in [24] that the boundary conditions are completely specified by the subgroup, and boundaries twisted by different cocycles for given subgroup can be adiabatically connected. Our results are further evidence in support of the observation in [24]. Figure 9: An illustration of the action of A v located at the physical boundary before and after they are shifted by a common group element g.

Conclusion
In this paper, we propose a slightly modified version of the calculation of entanglement entropy on lattice gauge theories, based on an observation locating a block diagonal structure of the reduced density matrix. This method explicitly reproduces the known results of the entanglement entropy of lattice gauge theories on closed surfaces.
We then apply the method to compute entanglement entropy when there are nontrivial gapped PBs. We show that the entanglement entropy has very subtle dependence on the group structure. As explained, the physics is intimately related to the interplay of condensed anyons at different boundaries. In the case of Abelian theories, the result has a simple closed form, but not so much so for more generic non-Abelian theories. It would be interesting to explore the physical implications of these results, particularly their connection to the structure of Frobenius algebra that underlies anyon condensation, and also modular invariants that they correspond to. We note that the analysis discussed in the current paper is applicable also for the TQD. The U (1) 3 cocycles would cancel in a very similar manner observed in section 4.8 where boundary 2-cocycles were canceled. It would be important to explore more general topological orders in this light beyond those describable by Dijkgraaf-Witten models. We note that some results on entanglement entropy in the presence of topological defects/interfaces have been explored in the context of Chern-Simons theories in the literature [33,34,35]. Topological boundaries can be considered as special topological defects/interfaces. How our results are related to existing observations should be explored in greater depth in a future publication.
Our next step would be to push these computations to higher dimensions. It is yet an open problem to have a complete classification of all possible topological boundary conditions in topological theories above 2+1 dimensions, let alone a unifying physical picture of these boundaries tantamount to the picture of anyon condensation. An understanding from the point of view of entanglement entropy should give new insights to solving the problem.
We also note that entanglement entropy can be used as a probe of higher form symmetries [36]. Our results might find fresh physical interpretations in terms of a connection to anomalies -or absence thereof -at the physical boundary.
We will leave these interesting and exciting questions to a forth-coming publication.
Note: We were notified of [37], which appeared in January when our paper were in preparation then, after posting our paper. In [37], it also explores the effect on the entanglement entropy in the presence of physical boundaries. In fact various situations considered in section 4 in the current paper have also been considered there, albeit following a different set of perspectives. In section 4.3-5 we considered the most general situations where the boundary conditions on the two (or more) boundaries of a manifold to be different. To our knowledge, this is the first instance it is considered in the literature. There are also some related results computing maximal/minimal entangled state on a cylinder in [38] mainly considering the Z 2 toric code. We hope that our methods would supply a useful alternative to the literuature.

A Open Ribbon operators and entanglement
Ribbon operators are defined in figure 10.
Let the ribbon operator F h,g act on the ground state. We choose a connected region containing one end of the ribbon as region R and the complementary region asR as shown in figure 11. We then compute the entanglement entropy between R andR.
We would like to compute the entanglement in the presence of these open ribbon operators, with one end lying in region R, and the other inR. We note that this has been considered before in [39], although the perspective and method adopted here is simpler. After the action of the ribbon operator F h,g , a fixed EB configuration A v i (g i ) does not correspond to a direct product state. Suppose on the two ends of the ribbon the A v operators are A v 1 (g 1 ) and A v 2 (g 2 ). Then after the action of F h,g g 1 , g 2 must satisfy g 1 g −1 2 = g. So we have g 2 = g −1 g 1 . We can use g 1 to label the state of region R and use region R F h,g Figure 11: Acting the ribbon operator F h,g on the ground state.
g 2 to labelR. Then for a fixed EB configuration, the state of the system can be written as |ψ = Here C stands for the conjugacy class, R stands for the irreducible representation of a representative element r C in the conjugacy class C, and N C is the centralizer of r C . u and v are two sets of integer parameters: (i, i ) and (j, j ). {c i } are the elements of the conjugacy class C. For each c i we fix a q i such that c i = q i r C q −1 i . Γ jj R is the matrix element of representation R.
Suppose the A v operator on the end of the ribbon is A v 1 (m). The the A v operator on the other end is A v 2 (q i n −1 q −1 i m). Using m and q i n −1 q −1 i m to label the state of R andR as previous, for a fixed EB configuration, we have We can observe that (A.5) Because of the orthogonality relation in the group representation theory, the above fomulation is a Schmidt decomposition. And we can still divide the EB configurations into G-orbits. So the entanglement entropy is where d i = |C| · dim(Γ) is the quantum dimension of the anyon type that corresponds to the ribbon operator F RC;uv . It matches with the result in [41].