On entanglement entropy in non-Abelian lattice gauge theory and 3D quantum gravity

Entanglement entropy is a valuable tool for characterizing the correlation structure of quantum field theories. When applied to gauge theories, subtleties arise which prevent the factorization of the Hilbert space underlying the notion of entanglement entropy. Borrowing techniques from extended topological field theories, we introduce a new definition of entanglement entropy for both Abelian and non-Abelian gauge theories. Being based on the notion of excitations, it provides a completely relational way of defining regions. Therefore, it naturally applies to background independent theories, e.g. gravity, by circumventing the difficulty of specifying the position of the entangling surface. We relate our construction to earlier proposals and argue that it brings these closer to each other. In particular, it yields the non-Abelian analogue of the"magnetic centre choice", as obtained through an extended-Hilbert-space method, but applied to the recently introduced fusion basis for 3D lattice gauge theories. We point out that the different definitions of entanglement theory can be related to a choice of (squeezed) vacuum state.


I. INTRODUCTION
Entanglement entropy has become an important tool for characterizing the correlation structure of quantum field theories [1][2][3], in particular with regard to correlations in space. In the latter case one presupposes that field degrees of freedom can be localized.
Gauge theories, however, feature a form of non-locality that prevents the strict localization of the so-called physical, as opposed to gauge-variant, degrees of freedom. For instance, in Yang Mills theories, including electromagnetism, the presence of Gauß constraints implies that one can compute the total electric charge contained in a region solely in terms of the electric flux across the region's boundary: no information about the bulk fields is needed.
Quantum mechanically, this non-locality is reflected in the fact that the Hilbert space of gauge-invariant states does not factorize into the tensor product of Hilbert spaces associated to a spacetime region A and its complement B. More precisely the algebra of gauge invariant observables does not factorize into the product of two commuting subalgebras each containing only operators supported in either A or B. Consequently, the definition of the entanglement entropy between one region and its complement requires further discussion, especially in the light of the privileged role gauge theories play in nature.
In quantum gravity this problem appears even more cogent, see also the recent discussions [4]. This happens not so much because a complete quantum theory of gravity is yet to be defined and agreed upon (in three dimensions one can actually argue for the opposite), but rather because of the very defining property of gravity: background independence. Indeed, because of background independence, which implies diffeomorphism invariance, the localization of regions and their separation into distinguished subsystems, when performed from within the theory itself, is already a thorny subject. To address the definition of entanglement entropy in a background-independent fashion, we advance a proposal which we believe sheds light also onto some of the issues encountered already within the standard gaugetheoretical framework.
To start with let us discuss the extant proposals for lattice gauge theories. In this context, essentially two approaches have been proposed for how to define entanglement entropy between two regions. Our work will further relate these two approaches and provide a new point of view.
The first approach [5][6][7][8], in particular put forward by Donnelly, is based on the embedding of the Hilbert space of gauge invariant states-which displays the non-local features discussed above-into an extended Hilbert space where gauge-invariance violations are allowed at the interface. This extended Hilbert space does factorize, allowing to define adjusted to the choice of vacuum. For states which have almost everywhere vanishing electrical field we should choose electric boundary conditions. Considering states which have almost everywhere vanishing magnetic flux-or in the gravitational language have almost everywhere vanishing curvature-we should choose instead the magnetic centre. These choices are designed to give finite results in the continuum limit, whose result is by construction already fully captured at the level of a finite (fine enough) lattice.
Furthermore, our definition of excitations is rooted in the analysis of the properties that regions with boundaries manifest under gluing. In fact, we will also argue that the process of cutting a system in two and the related definition of an extended Hilbert space procedure should be understood as dual to gluing. This brings into play techniques of extended topological field theory, i.e. topological field theory for manifolds with boundaries.
In particular, the choice of the BF vacuum will naturally lead us to consider a different basis for the Hilbert space, the 'fusion basis'. Importantly, thanks to a precise definition of its excitations, we can operationally specify regions by their excitation content. In this way we overcome the problems arising from defining a region independently of its content (something that would be in stark contrast with background independence), and we do so in a way that makes such a definition independent of the chosen regularizing lattice, thus directly avoiding the need of specifying the way we split it.
Indeed, our proposal comes to full fruition in (2 + 1)-dimensional gravity. This can be formulated as a gauge theory and moreover does not feature propagating degrees of freedom. Indeed it can be formulated as a topological theory, known as BF theory (from which the BF vacua are derived). This theory can be coupled to point particles which carry exactly the type of excitations considered in this paper. Therefore, regions will be specified by the point particles they contain. In this sense, our notion of entanglement theory characterizes the correlations between the excitations the regions contain.
We will focus our discussion to (2 + 1)-dimensional lattice gauge theories, and restrict for simplicity to finite gauge groups. A strategy to generalize the fusion basis, employed here, to (3 + 1) dimensions, can be found in [20]. The generalization to the Lie-group case requires the definition of the so-called Drinfel'd double for such groups, which can be found in [21,22].
The paper is organized as follows. In section II, we present the gluing of spatial manifolds with boundary. We then define the splitting procedure dual to this gluing as well as the notion of extended Hilbert space. In section III, we present the extended Hilbert space method of computing entanglement entropy and emphasize the relation with the observable algebra approach. In section IV, we introduce the set-up for lattice gauge theory and perform a change of point of view by focusing on the excitation content of the theory instead of the underlying lattice. This leads to the construction of the fusion basis, which will be key to our definition of extended Hilbert space. A new notion of entanglement entropy for lattice gauge theories is introduced in section V where several explicit calculations are presented. Finally, in section VI, we discuss the implications of this new definition for the case of 3D quantum gravity.

II. SURFACE GLUINGS AND SPLITTINGS, AND EXTENDED HILBERT SPACES
To define the entanglement entropy of a subregion, one first needs to specify how to associate the theory's degrees of freedom to it. In gauge theories, having to consider non-local gauge-invariant degrees of freedom, this leads to ambiguities.
The gauge-invariance condition manifest itself in terms of constraints, i.e. quantum versions of the elliptic equations a state must satisfy to represent valid initial data. In lattice gauge theories, this leads to Gauß constraints defined at the lattice nodes. The Gauß constraint at the node n involves all the links adjacent to it, links which carry gauge co-variant degrees of freedoms expressed in terms of parallel transports (holonomies) along open paths. As a consequence, a gauge-invariant wave function necessarily correlates the degrees of freedom across the links.
In turn, this prevents the splitting of the Hilbert space of gauge-invariant functions H into a tensor product H A ⊗H B , with the two factors associated to two complementary regions. The so-called extended Hilbert space procedure [7,8] avoids this task, by considering an extended Hilbert space H ext , in which the Gauß constraints are relaxed along the boundary interface between the two regions. More precisely, this interface is defined to be transversal to the links of the lattice, and a two-valent node is introduced on each link cut by it. The Gauß constraints are then relaxed for these two-valent nodes only. This defines an extended Hilbert space, which does factorize H ext = H A ⊗ H B in a straightforward manner.
In this work we will show that the extended Hilbert space procedure can be generalized using a different setup and also different sets of constraints. Furthermore this generalized procedure is deeply connected to the theory of extended topological field theories. There, indeed, one considers topological field theories on manifolds with boundaries 3 together with a procedure for gluing them to one-another. Splitting a manifold into two components thus arises as an inverse procedure.
In the following we will describe the main idea and start with the notion of gluing states, defined on spatial manifolds with boundary. Dual to this gluing one can define a splitting procedure and the notion of extended Hilbert spaces.

A. Gluing
To be concrete, we consider a theory where the (gauge co-variant) degrees of freedom are associated to the links of a graph Γ, as is the case in lattice gauge theories, where one has group elements g l associated to the links l of Γ.
Let Γ A and Γ B be two graphs embedded into the hypersurfaces, Σ A and Σ B respectively . We assume both Σ A and Σ B have boundaries and that the embedded graphs end at the boundaries by one or several open edges. The (gauge-covariant) wave functions defined on Γ A and Γ B respectively live in the kinematical Hilbert spaces H kin A and H kin B . Let then {C } A be a set of constraints which we require to be quasi-local, i.e. local e.g. with respect to the graph's nodes and faces and their adjacent structures, and H A the subspace of wave functions in the Hilbert space H kin A which satisfies these constraints. 4 Likewise for {C } B . To every constraint C , we assign a projector P C which projects onto the subspace of wave functions ψ satisfying the constraint: C ψ = 0.
As an example, we mentioned the Gauß constraint above, which acts at the nodes of the graph imposing gauge invariance. We consider however only the Gauß constraint acting at the internal nodes, that is everywhere except at the (one-valent) nodes on the boundaries 5 of Σ A or Σ B .
At the level of the surfaces, the gluing between Σ A and Σ B is obtained by identifying a portion of their boundaries. We denote the result of this operation Σ A∪B . At the level of the embedded graphs, it is analogously defined by connecting the links along which the gluing is performed. We denote the result Γ A∪B . Here we assume that the links ending at the two boundaries match under the gluing procedure. 6 Let us now define the gluing operation for two wave-functions and denote it with a star : H A ⊗ H B → H A∪B , where the wave functions in H A∪B have to satisfy the set of constraints {C } A∪B , which we specify presently. The gluing is defined in two steps.
Given two wave functions ψ A ∈ H A and φ B ∈ H B , consider first the (usual C-)product of wave functions ψ A · φ B , defined on the glued graph Γ A∪B . In general, this product wave function will not satisfy all the constraints {C } A∪B , which will include {C } A and {C } B but also further constraints that result from the presence of new internal nodes and faces in Γ A∪B . Nevertheless, the set of wave functions of the form ψ A · φ B will span the extended Hilbert space H A ⊗ H B =: H ext , and hence the subspace H A∪B can be identified with the set of subspace of wave functions which satisfy all the constraints {C } A∪B . Denoting the corresponding projector by P A∪B we finally define the star product as In the following we will denote for brevity, C = A ∪ B.

B. Splitting
Splitting is the 'inverse' operation of gluing. Given a surface Σ C and embedded graph Γ C we firstly have to introduce a boundary that splits Σ C into Σ A and Σ B and Γ C into Γ A and Γ B so that the gluing gives back the corresponding structures. (One might want to choose certain restrictions on which kinds of boundaries and graphs are allowed.) To define the splitting of a wave function in H C we are looking for an isometric embedding map 2) 3 Here we are working in a Hamiltonian framework, therefore 'boundaries' have to be understood as codimension 2 surface, which-in a covariant context-are usually called 'corners'. 4 Here we make the simplifying assumptions that the set of solutions to the constraints can inherit the inner product of H kin A . This happens if zero is in the discrete spectrum of the constraints. If this is not the case a new inner product needs to be constructed, see e.g. [23][24][25][26], and the following procedure needs to be amended accordingly. 5 Imposing a Gauß constraint for these one-valent nodes would mean to trivialize the dependence of the wave function on the group element associated to the link ending at the boundary. 6 This can be ensured by introducing marked points on the boundaries where the links are allowed to end. such that • E = id on H C . Note that this latter condition does not specify E uniquely, but we can demand that E maps H C into H C understood as a subspace of H ext . (Remember that H C can be identified with the set of wave functions in H ext satisfying all the constraints {C } C .) The embedding E, therefore, maps wave functions in H C , which does not allow a straightforward splitting, into an extended Hilbert space H ext H A ⊗ H B , which comes with a natural tensor factorization associated to the splitting of Γ C into Γ A and Γ B . Hence, to integrate out the degrees of freedom associated e.g. to Γ B , one first uses the embedding, and then traces over H B .

C. Flatness constraint
We have already mentioned the Gauß constraint of gauge theories. However, for our proposal the introduction of another constraint will be relevant. This is the flatness constraint, which acts at the faces of the graph and demands that these have a trivial holonomy, 7 i.e. a trivial magnetic flux through them. Again, we will have closed faces, necessarily internal to Σ, and open faces as well, necessarily including boundary components. We will demand the flatness constraints to hold for closed faces only.
At this point the reader might wonder why we are interested in the flatness constraints. Firstly, we can allow for curvature by introducing punctures, that is by removing disks from Σ and thus introducing boundaries around which the flatness constraint does not need to hold. Introducing sufficiently many punctures we can regain all curvature degrees of freedom. Therefore, our procedure is not over-restrictive. On the other hand, the introduction of flatness constraints allows to achieve a certain independence from the lattice by putting the focus rather on the punctures themselves, which provide the support for the excitations. Secondly, with regard to the gluing and extension process, the flatness constraints allow us to trade all Gauß constraint violations, appearing in the extended Hilbert space as defined in [7,8], for one flatness and Gauß constraint violation. The reason is that with the flatness constraints holding almost everywhere (except at the punctures) we can change the graph and its embedding without changing the physical content of the wave functions. Hence, in this way we can change the number of links crossing the boundary to just one link. In other words, the many local reference frames defined by the cut links are replaced by a global reference frame together with demanding a locally flat connection near the boundary. We will later see that this allows as to define the (generalization of the) 'magnetic centre' choice [9] in terms of an extended Hilbert space procedure.

D. Gluing of cylinders
Let us mention a further interesting application of the gluing procedure, which is crucial for the definition of the fusion basis and the notion of excitations. Consider a two-punctured sphere, denoted S 2 . It is topologically equivalent to a (hollow) cylinder. Let us denote by H S2 the Hilbert space of wave functions on a graph embedded in S 2 , satisfying the flatness and Gauß constraints for closed faces and inner nodes, respectively. Notice, that two copies of S 2 can be glued together giving another S 2 . Therefore, the state-gluing operation defines a multiplication map in S 2 : This equips H S2 with the structure of an algebra. More precisely, it turns out [27] that the multiplication provides a representation of the multiplication of the Drinfel'd double algebra D(G) of the gauge group G [28]. The (irreducible) representations of this Drinfel'd double algebra will play a central role in our proposal, in particular in the explicit construction of the fusion basis (cf. section IV F 2).

III. ENTANGLEMENT ENTROPY
We reviewed the issues arising when attempting the splitting of a Hilbert space H C of wave functions satisfying a set of constraints {C } C into a tensor product. Such a splitting can be performed by embedding the states in an extended Hilbert space H ext H A ⊗ H B for which some constraints are relaxed. This is described by an embedding map E : With a choice of embedding map at hand, we can define a notion of entanglement entropy for states in H C . To do so, we first use the map E to embed a given state ψ ∈ H C into H ext H A ⊗ H B , hence we define the reduced density matrix from which the entanglement entropy can be readily evaluated Notice that both D A and S A implicitly depend on E.
In [5,7,8], a definition of entanglement entropy was proposed for both Abelian and non-Abelian gauge theories by Donnelly. His procedure was the type we just described-often referred to as the 'extended Hilbert space' methodand made implicit use of a specific embedding map. In [9], CHR pointed out that (at least in the Abelian case) Donnelly's procedure agrees with their 'electric centre' prescription, but it was just one among other choices.
Here, we want to emphasize that, by choosing embedding maps different from Donnelly's, the extended Hilbert space construction can be generalized and is therefore not unique. In particular, the alternative procedure proposed here does reproduce CHR's 'magnetic centre' prescription. This at least for Abelian gauge theories, since we will see that the non-Abelian case necessarily includes also an electric component. Hence, in so doing, we provide a tighter connection between CHR's algebraic constructions and the extended Hilbert space procedure. Moreover, by explicitly providing an extended Hilbert space procedure matching the 'magnetic' centre prescription, we correct claims about its impossibility which have appeared in the literature [29].
In the rest of this section, we will describe in detail the contributions to the entanglement entropy, as defined by the extended Hilbert space procedure, along the lines of Donnelly [8]. While his analysis was based on a specific embedding procedure (corresponding to a choice of spin network basis for the Hilbert spaces involved), we will instead allow for generic embeddings and associated choices of basis. With this toolbox at hand, we will relate the extended Hilbert space procedure to CHR's observable-algebra-based definition [9]. It is left to the forthcoming sections, the task of introducing the details of the fusion basis for (2 + 1) dimensional lattice gauge theories [27], needed to parallel the magnetic centre choice, and the study of the corresponding embedding procedure and entanglement entropy.

A. Entanglement entropy from extended Hilbert spaces
Both the spin network basis and the fusion basis are indexed by representation labels, which we here will denote generically by ρ. Notice, however, that these are representations for different algebraic structures, namely for the group G and its Drinfel'd double D(G), respectively (cf. section IV F 2). As usual for basis-state labels, the ρ's encode the eigenvalues of a maximal set of commuting observables on the Hilbert space H C . This hints already at the connection to CHR's observable-algebra-based procedure as well as to more general choices of maximal sets of commuting observables.
We split the representation labels into three sets: {ρ A }, associated to region A; {ρ B }, associated to region B; and {ρ ∂ } associated to the boundary ∂A = ∂B. In the case of the fusion basis we will just have one ρ ∂ associated to the boundary. Also, in case we allow torsion excitation at the punctures-that is violations of gauge invariance there-we have representation space indices I associated to these punctures. These can be associated either to region A or B and we will therefore subsume them into the set of representation indices {ρ A } and {ρ B } respectively.
In the extended Hilbert space H ext = H A ⊗ H B we will have a basis that includes a doubling of the {ρ ∂ } labels to {ρ ∂A } and {ρ ∂B }. Furthermore, for each of these label sets, we have associated sets of representation-space labels 8 {I ∂A } and {I ∂B }.
Denote |ρ A , ρ B , ρ ∂ and |ρ A , ρ ∂A , I ∂A ⊗ |ρ B , ρ ∂B , I ∂B elements of an orthonormal basis of H C and H A ⊗ H B , respectively. The embedding map E is then given by The product is over the boundary elements. In the case of the fusion basis we will have only one boundary element and index ρ ∂ . In the case of the spin network basis any edge cut by the boundary is a boundary element.) Given a state the corresponding density matrix D A ψ , defined in (3.1), has the following structure: We see that D A is block diagonal, with each block labeled by a boundary-representation vector |ρ ∂A = ρ ∂ , I ∂A = I ∂ , and weighted by the probability distribution This distribution is constant over the I ∂ as a consequence of gauge invariance. On the other hand, the density matrix associated to each block (independent of I ∂ ) is Now, given this decomposition, one finds that the entanglement entropy (3.2) has three contributions [7] where • stands for the expectation value with respect to the classical probability distribution P (ρ ∂ ), and H(P (ρ ∂ )) for its Shannon entropy, B. Relation to observable-algebra-based entanglement entropy We now comment on the relation between this approach and the definition of entanglement entropy via the splitting of the observable algebra. CHR's original proposal [9] concerned only Abelian gauge theories. We will comment below on the non-Abelain generalizations.
Given the algebra of observables, O, associated to the gauge-invariant Hilbert space H C , one chooses a commuting subalgebra of observable, Z, associated to the boundary δA = δB. This commuting subset Z of observables serves as centre of a new, reduced, observable algebra O red , which is obtained by removing all the observables, which do not commute with the designated centre Z.
The choice of Z must be done in such a way that O red admits a splitting O red = O A ∪ O B into two mutually commuting subalgebras, which can be associated to the regions A and B, respectively. These subalgebras clearly have a non-vanishing intersection given by the centre, Now, H C (usually) provides an irreducible representation of the observable algebra O. By removing observables from O, one finds that the reduced algebra O red features superselection sectors on H C . These superselection sectors are precisely labeled by the eigenvalues {λ} Z for the observables in Z. This is because the original Hilbert space has by construction the structure The observables in Z can be also understood as boundary conditions, characterizing each of the superselection sectors. This interpretation physically explains the 'classical' behaviour of these observables noticed by CHR. See also the discussion in [30].
At the beginning of this section, we introduced the basis |ρ A , ρ B , ρ ∂ for H C . This basis immediately suggests one to choose the centre Z to be generated by the projectors P ρ ∂ onto the subspaces spanned by the |ρ A , ρ B , ρ ∂ with varying ρ A , ρ B but fixed ρ ∂ . This is equivalent to requiring the eigenvalues {λ} Z to be directly determined by the labels ρ ∂ ≡ {ρ ∂ } Z . Henceforth, with this choice in mind, we will replace the superindex {λ} by ρ ∂ .
The definition of entanglement entropy via the specification of a centre by CHR [9] did originally concern only the Abelian case. It can also be generalized to the non-Abelian case, albeit in two different manners. One choice corresponds to staying withing the algebraic framework based on gauge-invariant observables alone [10]. In this case one forms density matrices with a superselection structure as given by (3.10) (3.11) for each sector separately. The entanglement entropy for the entire system is then defined as where again • stands for the averaging with respect to P (ρ ∂ ).
We see that this result does not completely reproduce the extended Hilbert space procedure (3.8), as in (3.12) we do not have the term ln dim ρ ∂ appearing (this term trivially vanishes for the Abelian case). The source of the discrepancy is the following. In the extended Hilbert space procedure, the density matrices resulting from the embedding map E, have also a superselection structure. But this superselection structure is more refined: additional subsectors appear which are related to the internal indices of the representation I ∂ := I ∂A = I ∂B . Consequently, the density matrices resulting form the embedding procedure are effectively characterized by the following probability distribution and block density matrices (3.14) The entanglement entropy for density matrices with such a superselection structure is given by where here • P (ρ ∂ ,I ∂ ) denotes averaging with respect to P (ρ ∂ , I ∂ ) and • P (ρ ∂ ) averaging with respect to P (ρ ∂ ). This definition reproduces the splitting into three terms as in (3.8). This second choice of superselection structure, which includes the magnetic indices I ∂ , is not tied to the initial gauge invariant observable algebra, as the magnetic indices only arise after cutting the manifold. One can, however, argue that cutting the manifold one introduces a boundary, and that the magnetic indices should be part of the boundary data together with the ρ ∂ , characterizing (sectors of) wave functions defined on manifolds with boundary. In other words, one could argue that the splitting of a system into subsystem requires the introduction of additional information about the reference frames at the boundary, as encoded in the magnetic indices, which is needed to perform a consistent gluing.
Thus the first two contributions to the entanglement entropy in (3.15) are resulting from the superselection structure and are thus due to the classical probability distribution (3.13). Indeed it was conjectured by CHR and proven by [10,31] that only the third contribtion in (3.15) gives the so-called distillable entropy, which is defined to be the maximum number of Bell pairs that can be extracted by a so-called entanglement distillation. The latter process involves a choice of (local) operator algebra, which in [10,31] is based on the reduced operator algebra O red . Thus the notion of distillable entropy also depends on the choice of reduced operator algebra or alternatively boundary conditions.
From a physical standpoint, what all this discussion is reminding us is that the concept of entropy is coarse-graining, i.e. observer, dependent. By varying the amount of information we know, or conversely, we would like to know about a system, we calculate different entropies. This can be summarized in the statement, that entropy is an epistemological quantity. And in sophisticated enough situations, also the entanglement entropy is such.
As mentioned, the extended Hilbert space procedure was first proposed using spin network functions [5][6][7][8]. Here the representation labels ρ ∂ characterize the eigenvalues of electric flux operators associated to the links that are cut by the boundary (for non-Abelian gauge theories one can take Casimir operators formed from the electric fluxes associated to each such link, see [10]). Thus, in the corresponding algebraic definition the centre is formed by these electric operators.
In this paper, we focus on the extended Hilbert space procedure for the fusion basis, where ρ ∂ characterizes a so-called closed-ribbon operator along the boundary between the two regions A and B. In the case of an Abelian theory (and considering only gauge-invariant wave functions), this ribbon operator reduces to a Wilson loop. Thus, in this case, the centre is given by a 'magnetic' operator. For non-Abelian gauge theories, however, the closed-ribbon operator measures also an electric excitation, related to the total flux of the electric field flowing out of the enclosed region. Note that this can be non-trivial even for completely gauge invariant wave functions. Thus, for non-Abelian theories, the magnetic centre gets naturally enlarged by a further 'electric' operator.
Finally, we now turn to the introduction of the fusion basis.

IV. FUSION BASIS FOR LATTICE GAUGE THEORIES
In this section, we review the construction of the fusion basis for a (2 + 1)-dimensional lattice gauge system. For a more extensive treatment we refer the reader to [27]. For simplicity, we will assume that the gauge group is a finite group G. We will also fix the topology of the underlying two-dimensional hypersurface to be spherical (S) with possibly p punctures present, i.e. Σ S p .

A. Hilbert space HΓ
Let Γ be a graph embedded in S. To start with, assume that the graph has no open links, i.e. no links ending at one-valent nodes. The graph gauge-connection is defined by associating a group element to every (oriented) link of the graph corresponding to the holonomy along the link. The Hilbert space H Γ is spanned by the functionals ψ : G L → C on the space of holonomies, where L denotes the number of links in Γ. The Hilbert space H Γ is equipped with an inner product defined as Gauge transformations are parametrized by {u n } n ∈ G N , where n denotes a node of Γ and N the number of such nodes. A gauge transformation acts on a holonomy configuration {g l } ∈ G L as where s(l) and t(l) denote the source and the target nodes of the link l, respectively. Gauge invariant functions are functions invariant under this gauge action. This defines a subspace H G Γ ⊂ H Γ of gauge invariant functions in the Hilbert space H Γ . H G Γ inherits the inner product (4.1) from H Γ . The gauge-invariance condition is encoded in the following Gauß constraints (or projectors), associated to the nodes n: where s and t index the links for which n is a source and target node, respectively, while l indexes the remaining links.

B. Basic operators
Considering the configuration space to be the space of group holonomies, we have available two kinds of operators, namely Wilson loop operators and translation operators.
Wilson loop or (closed) holonomy operators W f γ act as multiplication operators on states ψ({g}). Given a function f : G → C and a path γ which coincides with some oriented and connected path along the links of Γ, the action of W f γ is given by where h γ = g ln · · · g l1 for γ = l n • · · · • l 1 . Note that for W f γ to commute with the Gauß constraints, γ has to be a closed path and f a class function.
Translation operators T k [H] act by finite translations and therefore correspond to an exponentiated version of flux operators. We can define both a left and a right action for such operators. We choose to work with left translation operators whose explicit action is given by (4.5) The action of T k [H] typically induces violations of the Gauß constraint at the target node of l k , i.e. t(l k ). However, the node at which the Gauß constraint violation occurs can be moved at will. For this, one can parallel transport the to-be-translated argument g k from its target node t(l k ) to some other node n along a path 9 γ, apply the translation in the frame of n, and then transport the resulting holonomy back. We denote these operators T k,γ [H] and their action reads where h γ was defined above as the holonomy along the path γ. The violation of the Gauß constraint induced by the translation now appears at the node n. This property will be important in the forthcoming construction where we will combine these basic operators in order to obtain so-called 'ribbon operators'.

C. Shift of viewpoint
Let Γ be a graph embedded into S. Γ being planar, we can unambiguously identify its plaquettes or faces. The shift of point of view we propose relies on the assumption that Γ carries excitations located at the faces. As we will explain presently, these excitations have to be understood with respect to a given vacuum.
First, we consider curvature excitations since they are naturally carried by the faces of Γ. Indeed, they are characterized by the amount of curvature carried by every face, defined as the trace of the holonomy surrounding the face. In the electromagnetic case these are precisely the magnetic fluxes.
Then, we consider torsion excitations, that is violations of the Gauß constraints (4.3). In the electromagnetic case these excitations correspond to the presence of non-vanishing electric charges. Being associated with a Gauß constraint violation, these excitations are a priori located at the nodes of Γ, and not at its faces as we desired. To obviate this problem, we introduce extra links and nodes. More precisely, we introduce one new link and one new node for each face (in the context of combinatorial quantization of Chern Simons theory, this structure is called a 'cilium'). For a given face this new link starts at some node on its boundary and ends at a new one-valent node placed in its interior. Henceforth, we refer to these nodes as 'end nodes' {n e }, and to all others nodes as 'internal nodes' {n i }. The valency of an internal node is strictly bigger than one. In the same spirit, we call the links adjacent to the end nodes 'open links'. Figure 1 depicts such a construction in the case of a lattice with square faces.
The result of this construction is an extended graph Γ which leads to a new Hilbert space H Γ equipped with an inner product of the same form as the previous one, see (4.1).
As it will become clear later on, allowing for torsion excitations is a necessity in the case of non-Abelian gauge theories, even in the case we do not allow them at the lattice scale. We restrict our focus on the subspace of H Γ constituted by wave functions which are gauge invariant at all internal nodes, but not at the end nodes. This defines a new Hilbert space, H p , where p stands for the number of end-nodes in Γ, which by construction coincides with the number of its faces, too: Note that the Hilbert space H p is unitarily equivalent to the subspace of wave functions in H Γ which are gauge invariant at all nodes where one does not attach an open link. In other words, we can map the torsion excitations from H p to H Γ by associating them with the nodes to which one attaches an open link. For the example of the lattice depicted in Figure 1 gauge invariance violations at almost all nodes can be taken into account in the Hilbert space H p . Furthermore one can also generalize the definition of H p , allowing more than one open link to end in a given face [27]. This allows to take into account all possible gauge invariance violations, starting from an arbitrary graph Γ. The change of point of view we adopt here can be made more explicit by placing a puncture in the middle of each face. More precisely, instead of thinking of a lattice embedded in S 2 and allow for some excitations, we can directly imagine a graph embedded onto a punctured sphere. The punctures act as defects which are the only possible support for both curvature and torsion. In this case, we can map a lattice with p faces to a graph embedded on a p-punctured sphere. The face holonomy becomes the holonomy surrounding the puncture while the open edges now go from a node of the graph to one-valent node sitting at the puncture. This correspondence is detailed in the next paragraph.

D.
From graph to punctures, and the flat vacuum The notion of excitations is bound to a notion of vacuum, which here is the state without any curvature and torsion excitations. This vacuum is a gauge invariant state peaked on flat connections, also known as BF vacuum. 10 Now, imagine that one wishes to describe configurations of (continuum) connection fields that are everywhere flat, except at a pre-defined number of points. Including torsion excitations, we need to extend these points to so-called punctures, that is infinitesimal disks with a marked point on their boundaries. The connection degrees of freedom can now be encoded in a Hilbert space H p as described above. Note however that the precise choice of graph does not matter. For instance due to the local flatness of the connection, we can deform links of the graph, as long as we are not crossing over a puncture.
Moreover we can even allow for graph refinements, that is add links, so that we have additional faces, that do not contain a puncture. In this case we just need to make sure that all allowed wave functions prescribe vanishing curvature at each closed face (i.e. at each face with no associated puncture) and, similarly, that they are gaugeinvariant at all internal nodes. In conclusion, we are allowed to change of graph, as long as it is sufficiently fine to (i) capture the first fundamental group of the punctured sphere, and to (ii) allow at least one connected path between any pair of punctures. See [27] for the precise transformation maps.
Later, we will also introduce 'ribbon' operators whose action on the states does not depend on the particular choice of underlying graph either. These are the operators operators which will be used to characterize the fusion basis.
Let us emphasize that this is a useful viewpoint one can adopt which makes explicit the connection to topological field theory with defects. However, although convenient, it is not necessary, and one can also proceed by having the usual fixed lattice in mind.

E. Holonomy basis for Hp
We now construct a holonomy basis of the Hilbert space H p . As the name suggests, this basis is designed to diagonalize holonomy operators. These operators are demanded to be based on paths which start and finish at the end nodes. To define a maximal set of such holonomy operators, we join the two following subsets: i) G-holonomies-First, we single out one end-node and call it the 'root node', n r . We call the face enclosing this root node the 'outer face'. We then need to choose a set of paths from the root node to each of the other end nodes {n e }. For this we pick a (connected) spanning tree in Γ denoted T . Such a tree uniquely determines a path P n from the root node to any other node n, and a fortiori also to the end nodes of Γ . The set of G-holonomies {G n e } is defined as the oriented product of holonomies following the paths P n e : This set automatically fixes all holonomy between pairs of end nodes along paths supported on T .
ii) H-holonomies-The second set is constituted of holonomies {H n e } based on closed paths {L n e }, going anticlockwise along the boundary of every face containing a puncture n e (all the others being trivial, anyway) and starting at the end-node n e associated to the face itself, Note that in order to obtain a maximal set of holonomies, it is not necessary to include the one around the outer face as long as we include the holonomies around all the other faces. Thus, a basis wave-function turns out to be labeled by (p − 1) pairs (G n e , H n e ) ∈ G 2 . Denote it ψ {G n e ,H n e } . Figure  2 depicts an example of such a construction. The wave functions can finally be written in a fully covariant form as a product over delta functions prescribing the G and H-holonomies. For the sake of clarity, let us look at the minimal examples of a lattice with two faces. Replacing the faces by punctures, this corresponds to considering an embedded graph on the two-punctured sphere. Since the two-punctured sphere S 2 is topologically equivalent to a cylinder, we have the following gaphical correspondence: with the marked point at the bottom puncture chosen as the root node. Applying the previous prescriptions, the gauge covariant form for the holonomy basis states on the two-faces square lattice (or two-punctured sphere S 2 ) is given by where we have chosen a particular normalization that will turn out to be convenient later on. All delta-functions are Kronecker deltas, whose value is either zero or one.

F. Fusion basis and ribbon operators
The holonomy basis {ψ {G n e ,H n e } } diagonalizes holonomy operators that are not gauge invariant at the end nodes and for this reason it is for now quite involved to specify a complete and independent subset of fully gauge invariant wave functions. Therefore, we first aim to find a (maximal) set of gauge-invariant operators and hence the basis which diagonalizes it.
Starting from the holonomy basis, the previous remark suggests that we should include into the set of gauge invariant operators the conjugacy class of the holonomies {h n e }. Let us for instance consider a lattice and two faces associated with the end nodes n e 1 and n e 2 . We denote h n e 2 ∪n e 1 the holonomy surrounding these two faces. If we have a non-Abelian group G, knowing only the conjugacy classes C 1 and C 2 of the two holonomies h n e 1 and h n e 2 , will in general not determine the conjugacy class of h n e 2 ∪n e 1 uniquely. Therefore the conjugacy class of the holonomy going around two faces generally encodes more information than the one provided by the conjugacy classes of the individual faces. It turns out that, knowing the individual conjugacy classes, the set of conjugacy classes one can obtain for the holonomy around the two faces is determined by so-called fusion rules.
We propose to construct a basis which relies on the notion of fusion sketched above. For this reason, we refer to it as the fusion basis. The fusion basis diagonalizes a hierarchical set of (gauge invariant) operators, detecting the conjugacy classes of loop based holonomies. This hierarchical set is described by a so-called fusion tree. We choose it to be rooted and binary (i.e. with three-valent internal vertices) such that the end vertices of this tree are associated to the faces (or punctures) together with their corresponding end nodes, and the root of the tree is associated to the outer face with the root node n r . The combinatorial structure of the fusion tree determines which faces (or loops, or punctures) and in which order, are fused to form larger ones. Thus, the fusion tree determines for which hierarchical merging of loops one considers the associated closed holonomies. As an extra condition, we require the set of loops underlying the closed holonomies not to cross each other.
The hierarchical set of loops { } defined above prescribe gauge invariant functionals {f (g )} which detect the conjugacy classes and therefore capture the curvature (or magnetic) degrees of freedom. In particular, this defines Wilson loop operators {W f }. However, we would also like to have operators that characterize the torsion (or electric) degrees of freedom. Indeed, even if we consider completely gauge-invariant functionals without torsion degrees of freedom for the original faces, we might have 'emergent' torsion degrees of freedom which arise when applying the fusion scheme described above. This is the one reason why torsion excitations might appear under coarse graining [17,27,[32][33][34]. This feature is again characteristic of non-Abelian groups and such effective torsion charges have been named Cheshire charges in [32]. Conveniently, the torsion degrees of freedom can be captured with operators based on the same hierarchical set of loops, as the one used for the curvature degrees of freedom. The difference is that these operators include the action of translation operators.
By putting together these two kinds of operators, we obtain the so-called ribbon operators introduced (in a slightly different form) by Kitaev [35]. These operators measure both curvature and torsion excitations. We are now ready to review their construction.

Closed ribbon operators
Let us now introduce the closed ribbon operators which are diagonalized by the fusion basis. We will see later that we can also define open ribbon operators. In order to describe the action of the closed ribbon operators, we will focus on a simple example, the general case as well as more details can be found in [27].
Consider the piece of graph displayed in Figure 3. There is a loop based holonomy given by g = g 2 g 1 where g 2 is the composed holonomy which goes from n 2 to n 1 while g 1 goes from n 1 to n 2 . We are going to define the action of a directed closed ribbon operator along this loop . By convention, the ribbon operator is drawn to the right (with respect to the orientation of the ribbon) of the associated loop . For the purpose of describing the action of the ribbon operator we have to choose an initial node, which here will be n 1 i.e. the target node of the link carrying h 1 . However, the action of the closed ribbon will eventually not depend on this choice.
Firstly, consider the action of the following (auxiliary) operator, parametrized by (G, H) ∈ G × G: which is a combination of a holonomy operator-the wave function is multiplied by δ(G, g 2 g 1 )-and a series of translation operators acting on the links crossed by the ribbon. The translation parameter is always parallel transported along the loop to the node n 1 . Both the holonomy action and the translational action a priori violate gaugeinvariance at this node. Furthermore, as the ribbon contains a translational part, it might induce (or modify existing) curvature, at the face which includes the holonomy h −1 1 g 2 h 2 . Indeed, this holonomy undergoes the shift On the other hand, note that the holonomy h −1 2 g 1 h 1 stays invariant.
Therefore the auxiliary operator R[G, H] might induce a change in curvature, via the action on the holonomy h −1 1 g 2 h 2 , as well as violations of the gauge invariance at the node n 1 . However, a closed ribbon operator is expected to only measure excitations, not to induce them. To obviate this problem, on the one hand we demand that G commutes with H, that is we require H to be in the stabilizer group N C of any representatives of the conjugacy class C of G. This prevents any curvature modifications. On the other hand, to deal with the gauge-invariance violation, we group-average the resulting wave function over the gauge action at the node n 1 . This amounts to consider the averaging over the adjoint action of the group on the parameters (G, H) of the ribbon. Putting everything together, one can show that the group averaged operator only depends on the conjugacy class C of G and a conjugacy class D of H ∈ N C : where the normalization factor |N D | corresponds to the cardinality of the stabilizer group N D of H ∈ N C . As expected the action of the operator K[C, D] does not depend any more on the choice of auxiliary node (which in this example was n 1 ). One can further show that the precise positioning of the ribbon with respect to the graph Γ is immaterial and the only information that matters about its position is topolgoical-i.e. it only matters how the ribbon winds around the locations of the excitations, that is the punctures. We have thus defined the closed ribbon operator K[C, D]. It is already clear that the closed ribbon operator projects onto states peaked sharply on the conjugacy class C for the loop based holonomy g = g 2 g 1 . We would also like to achieve a projector property with respect to the parameter D, prescribing the translational action of the ribbon operators on holonomies transversal to the loop holonomy. Usually, the diagonalization of a translation operator requires some sort of Fourier transform. Indeed, we define where R denotes a unitary irreducible representation of the stabilizer group N C , χ R the corresponding character, and d R its dimension. We name these newly defined operators K[C, R] 'charge ribbon operators'. They are projectors: In summary, we have obtained closed ribbon operators K[C, R] which measure the excitation content of the region enclosed by the ribbon. This excitation content is characterized by two parameters: a conjugacy class C of G and an irreducible representation R of the stabilizer group N C . The conjugacy class C describes the curvature excitations, whereas R measures the torsion. In more concrete physical terms in the gravitational context, we can understand the punctures as point particles coupled to (2 + 1) gravity. In this case, C encodes the mass of the particles, while R the component of its spin (projected along the internal direction defined by the curvature). For Yang-Mills theories, on the other hand, C is a measure of magnetic flux, while R can be seen to measure the integrated flow of electric field into the region enclosed by the ribbon. This is because the translation operators are exponentiated versions of what would be an electric flux operator if G were a Lie group.

Fusion basis
A set of closed ribbon operators {K β [C, R]} is mutually commuting as long as the ribbons do not cross each other. The fusion basis diagonalizes exactly a certain choice of such mutually commuting closed ribbon operators. This leads to a hierarchical set of ribbons. Indeed, first we consider the set of ribbons around the basic faces (or punctures), excluding the root face. These define the basic excitations. One then fuses two excitations by considering ribbons around fused faces or punctures. In each step one fuses two excitations, which can be either basic ones or excitations resulting themselves from a fusion. One proceeds until there is only the outer face or root puncture left. Since we consider a sphere the ribbon around the root puncture agrees with the ribbon around the remaining punctures modulo orientation. Notice also that the faces do not have to be neighbouring with respect to a particular choice of underlying graph (this is why it is more powerful to get rid of the graph as we advocated above), but it is important that in the final set of closed ribbons, no two ribbons should cross each other. The choice of fusion scheme can be encoded in a fusion tree, where the end vertices of the tree correspond to the end nodes of the graph, that is to its faces or punctures, and the root of the tree correspnds the root of the graph, that is to its outer face. The trivalent vertices of the tree encode the fusion of two excitations into a new one. The edges of the fusion tree are labeled by pairs (C β , R β ) where β ∈ {1, · · · , 2p − 3}. These labels determine a set of fusion basis states, namely those fusion basis states the closed ribbon operators K β [C β , R β ] project onto.
Let us now construct explicitly such fusion basis states. We start with a holonomy basis state ψ G,H whose support is a two-punctured sphere S 2 or equivalently a graph with two faces. Here the basis state ψ G,H is labeled by the pair (G, H) ∈ G × G. We explained in section II, that the gluing procedure for wave functions on S 2 reproduces the multiplication map of a well-known algebraic structure, namely the Drinfel'd double D(G) of the group G. In particular the Drinfel'd double admits a basis {[G, H]} labeled by couples (G, H) ∈ G 2 (cf. appendix A for a brief review). One can then define elementary excitations by demanding that the corresponding wave functions are stable under such a gluing procedure [36]. This leads to the identification of elementary excitations with the irreducible representations of the Drinfel'd double D(G), whose construction we now briefly briefly summarize [37][38][39]. Following an inducedrepresentation type of construction, one can show that the irreducible representations ρ of D(G) are labeled by couples ρ = (C, R), where-as above-C is a conjugacy class of G and R an irreducible representation of the stabilizer group As S 2 is topologically equivalent to a cylinder we introduce the following graphical notation  We see, however, that these operators K[C, R] do not suffice to fully characterize the fusion basis, because of the presence of further basis label associated to the punctures (i.e. I and I ) . One can introduce projection operators P α [ρ, I ], whose purpose it is to project onto a fusion basis state carrying the label (ρ, I ) at the puncture α. Therefore, the closed ribbon operators, together with these projection operators, give a complete set of commuting operators, characterizing the fusion basis. With the irreducible representations at hand, we can make the notion of fusion of the basisc excitations more precise. As we have seen, the elementary excitations are described by irreducible representations ρ of the Drinfel'd double, so that the fusion of two excitations is described by its recoupling theory: For notational convenience, we will assume that N ρ3 ρ1ρ2 either vanishes or is equal to 1 (i.e. that the irreducible representations of D(G) form a multiplicity free fusion category), however the following derivations still hold without this assumption. With a choice of basis (and phases) for the representation spaces, the decomposition is described by Clebsch-Gordan coefficients satisfying and which can be graphically represented as Relaxing the multiplicity-free assumption would lead to an additional multiplicity index for the Clebsch-Gordan coefficients.
We have now all the ingredients to define the fusion basis for the general case. Let ψ {Gα,Hα}α be a holonomy basis state whose support is a p-punctured sphere. Since the number of punctures is p, the basis state is labeled by (p − 1) pairs (G α , H α ) ∈ G 2 . Each one of these pairs corresponds to a basis excitation and can be thought as labelling a cylinder state ψ Gα,Hα . To define the fusion basis state, we first need to perform the transformation to the [ρ , I I]-picture on each of these (p − 1) pairs respectively associated to (p − 1) punctures: Here the index I α is associated to the corresponding (non-root) puncture α, whereas the index set {I α } α are all associated to the root puncture. Thus the root puncture carries the tensor product over all representations ρ α . To make the fusion explicit, we decompose this tensor product by using a recoupling (or fusion) scheme encoded in a choice of fusion tree. That is, for every three-valent vertex of the fusion tree we apply a Clebsch-Gordan coefficient to the states (4.23). The contraction of all Clebsch-Gordan coefficients according to the fusion tree leads to a tensor C Equivalently, we can identify (4.23) with a product of states on the cylinder and write the fusion basis as This leads to the following encoding of the fusion basis in a graphical representation, which also includes the fusion tree: (4.26) Here the (p − 1) upper cylinders are associated to the (p − 1) punctures and are connected to each other via Clebsch-Gordan coefficients. This basis can be shown to be orthonormal and complete [27], and to diagonalizes the closed ribbon operators supported on the loops associated with the relevant fusion tree: Moreover, with the above graphical notation, it is clear that the usual lattice-based representation has been abandoned in favour of a representation relying exclusively on the excitations and the way they fuse together. Therefore it is natural at this point to define a region not so much in terms of the underlying lattice but only in terms of the excitations it contains.
Note that a different choice of fusion tree would lead to different fusion basis. The fusion basis given above is based on a particular such choice of tree. In the following, see also appendix B, other choices will turn out to be more relevant.

Fully gauge-invariant wave functions
Fully gauge-invariant wave functions can be also easily described in terms of the fusion basis. The gauge invariant projection at the puncture α, implies that it carries a trivial representation labels R α = 0. This trivial representation label R α = 0, however, still allows the index I α = (i, M ≡ 0) to range among the values of i labeling the elements of the quotient Q Cα = G/N Cα . The gauge-invariant projection induces, however, an averaging over the elements of Q Cα , labeled by i. Therefore, the gauge-invariant projection of a fusion-basis states is effectively labeled only by the conjugacy classes C α . Note, however, that in the non-Abelian case fully gauge-invariant basis states can actually be labeled by ρ β = (C β , R β ) which include non-trivial representations R β for β's associated to inner edges of the fusion tree.

Open ribbon operators
We have already introduced the closed ribbon operators, which measure the excitations without changing the excitation content.  The ribbon operator R[G, H] is parametrized by (G, H) ∈ G 2 and combines (as before) a holonomy operator part, acting on the holonomy parallel to the ribbons, and a translation operator part, acting on the holonomies crossed by the ribbon: The parallel transport for the translation part insures that gauge invariance is preserved at the internal nodes. Furthermore, the action is defined such that the curvature is changed only for the faces that include the source and target end-nodes of the ribbon.
Let us define a vacuum state as being a state without excitations, that is, the curvature is vanishing for all faces (or punctures). Vanishing torsion, on the other hand, means that the states are completely spread over the holonomies going from one puncture to another. Thus the (BF) vacuum state is given in the holonomy basis and up to normalizations by 29) where e denotes the unit element of G. Applying a ribbon operator to the vacuum state, we create curvature and torsion excitations at its ends. Because we can create such excitations only in pairs, they are of a quasi-local nature. We can also ask for charge ribbon operators that would create basic excitations, labeled by the Drinfeld double algebra representations ρ = (C, R). In fact, applying the generalized Fourier transform, we obtain the operators which indeed generate the fusion basis for the two-punctured sphere from the vacuum state From the relation between the fusion basis and the cylinder states given in equation (4.25) it should be evident that by appropriately acting with these charged ribbon operators on the vacuum, one can indeed generate the whole fusion basis. An important related fact is that it is possible to define a "lengthwise" product of open ribbon operatorsappropriately sharing a start and end-puncture-which leads to a new open ribbon operators supported on the combined path. This multiplication reflects the gluing operation and satisfies the same Drinfel'd double algebra (see [27] for details).

V. ENTANGLEMENT ENTROPY IN LATTICE GAUGE THEORIES
Having introduced the fusion basis, we now have all the ingredients to define a new notion of entanglement entropy through the procedure of section III. There, we assumed to have a basis labeled by representation indices {ρ}, and that the system under scrutiny was divided in two regions generically associated to a set of representation labels {ρ ∂ }. Now, in the case of the fusion basis, and assuming that there is only one connected boundary, this set of labels can always be reduced so that it includes only one representation ρ ∂ . This representation label describes the outcomes of the closed ribbon operator going along the boundary between regions. The extension procedure described in section II, introduces an extended Hilbert space H ext = H A ⊗ H B , which factorizes into two Hilbert spaces. These two Hilbert spaces correspond to the two systems one obtains after splitting the surface Σ along the boundary, see the discussion in section II.  5. For a given fusion tree we identify the region A as a set of punctures and B its complement. The splitting is performed by cutting a cylinder of the fusion tree. This cut requires the introduction of additional punctures. Note that the fusion tree employed here is the same as the one appearing in the alternative fusion states defined in Appendix B. Figure 5 represents the cut into two regions in the picture using punctures and the fusion tree. In the usual lattice picture the cut proceeds along the boundary of the plaquettes. More precisely we can imagine to double the Wilson loop around e.g. the region A into two loops which go closely parallel to each other. The two Wilson loops are connected with one 'small' link, whereas the area between the two loops carries flat connection. (Demanding flat connection for this area and gauge invariance at the additional nodes, we can uniquely map the state on the original graph to the graph with the doubled Wilson loop.) The cut proceeds then in-between the two Wilson loop and cuts the link connecting the two loops. Let us mention some key differences between this extension procedure and the one which uses spin network states as in [5,7,8]. The use of the fusion basis emphasizes the role of the excitations with respect to the BF vacuum, which in turn prescribes gauge-invariant flat connection. In (2 + 1) dimensions this means that only the position of the punctures matters, not the graph itself. This feature is particularly important for the application to (2+1) dimensional gravity, where the BF vacuum, with defect excitations describing particles, give the physical states [40][41][42].
The two procedures can be also compared in how 'big' the extension of the Hilbert spaces is. This can be quantified in terms of how many constraints are violated in the extended Hilbert space. In the case of the spin network basis this includes all the Gauß constraints at the two-valent nodes that result from the boundary cutting links. In the case of the fusion basis this includes only one Gauß constraint and one flatness constraint. In the picture described above the Gauß constraint needs not to hold anymore for the link which is cut into two. And the flatness constraint that is violated in the extended Hilbert space is the one between the two Wilson loops arising from doubling the Wilson loop along the boundary. This is, nonetheless, a sort of minimal choice-one that can, importantly, always be made. A more general splitting can be introduced also for the fusion basis. This corresponds to choosing somewhat arbitrary more marked points along the boundary of the two regions. In this way, more Gauß constraint violations will be added (but no extra curvature violation, see [27]). From the observable algebra perspective, this corresponds to declaring observable a series of gauge-variant holonomies along a set of paths which partition the boundary. This extension will only introduce additional vector-space indices (i.e. new indices 'next to' the (I, I ) in V ρ ∂ ). Later, we will briefly discuss what kind of consequence this has for the entanglement entropy.
We are now going to explicitly calculate the entanglement entropy for some simple choices of states and regions. To remind the reader, the entanglement entropy has three contributions where P (ρ ∂ ) is the classical probability distribution while • and H(P (ρ ∂ )) denotes the average with respect to P (ρ ∂ ) and its Shannon entropy respectively.

A. Fusion basis states (and BF vacuum)
We start with a fusion basis state on the p-punctured sphere S p . Its expansion in the fusion basis of course gives Here β labels the edges of the fusion tree and α its endpoints. We partition the punctures into two sets A and B.
We choose a fusion tree such that the A and B sets are only connected by one fusion tree edge labeled by ρ ∂ = ρ 0 ∂ . That is we are only considering fusion basis states from a basis characterized by a tree which is 'compatible' with the prescribed splitting. Such a basis always exists (and in fact, there are many).
In this case the classical probability distribution P (ρ ∂ ) is peaked on one particular value P (ρ ∂ ) = δ(ρ ∂ , ρ 0 ∂ ), thus the associated Shannon entropy vanishes. The density matrices D A (ρ ∂ ) are defined to vanish for ρ ∂ = ρ 0 ∂ , and the density matrix D A (ρ 0 ∂ ) has only one non-vanishing entry equal to 1 on the diagonal, and therefore give no contribution to the entanglement entropy (5.1). Therefore, we are left with the middle term in (5.1) Notice that we find a vanishing entropy for the BF vacuum state, as in this case dim ρ 0 ∂ = 1. This agrees with the result found in [9] for Abelian gauge theories with the magnetic centre choice on the BF vacuum state. For Abelian structure groups, ρ = (C, R) is labeled by a group element (as C = {g}) and an irreducible representation of its stabilizer, that is of the whole group. This irreducible representation is-the group being Abelian-one-dimensional. Hence, in this case dim ρ = (dim C)(dim V R ) = 1, and the entanglement entropy vanishes for (compatible) fusion basis states.
We can also consider gauge invariant projections of fusion basis states, as defined in section IV. We can consider these states both in the Hilbert space H p which allows for torsion excitations at the punctures, and its gauge invariant projection, where torsion excitations do not appear for the punctures (but can-in the non-Abelian case-appear for internal edges of the fusion tree). In both cases the result is the same as in (5.4). The only difference could have been in the contributions from the density matrices D A (ρ ∂ ), but these do describe pure states also after the action of gauge-averaging projectors. This is because such projectors act locally within one single region.
Let us shortly come to the extension mentioned at the end of this section's introduction. This extension consists in a generalization of the Hilbert space H p and of the corresponding fusion basis to the case where more than one link is allowed to end on a given puncture (in this case the puncture that is identified with the boundary between the two regions, see [27]). It turns out that the corresponding fusion basis is still labeled with the same representations, the only difference is in the vector space indices I at the extended puncture: for each additional marked point accompanied with a link ending at the puncture the index range is multiplied by |G|, the order of the group. Associated to this generalization of the fusion basis we can also consider a generalization of the gluing and extension procedure. Basically we can decide by how many graph links we wish to connect region A and region B. Note that this reintroduces a graph dependence 11 that previously we fixed by using an equivalence relation between states, that allowed us to change the underlying graphs. This enabled us to always reduce to the case that region A and region B are connected by only one link.
Using this generalized procedure the adjustment of the entropy formula is very simple. The result is that for each additional marked point the entanglement entropy increases by ln |G|. Interestingly, it turns out that the following relations between dimensions hold (see e.g. [27]) where the first is a dimension of the Drinfel'd double seen as a vector space spanned by the basis {[G, H]}, the second is the cardinality (order) of the group G, and the last one is again a dimension of a vector space, i.e. of V ρ . The last term in the above equality is also known as the (square of the) 'total quantum dimension' of the fusion category given by the irreducible representations of D(G). Due to the lack of fonts for a capital 'D', we will call this quantity, i.e. the total quantum dimension of D(G), simply Ω: Hence, we find that for a (compatible) extended fusion basis state, with m marked points 12 at the boundary puncture, itself labeled by ρ ∂ , the entanglement entropy amounts to Notice, that by using such a graph-dependent formula, one obtains also a non-vanishing contribution for the BF vacuum state. This can be understood by realizing that the extended Hilbert space allows now a refined information on the gauge connection along the boundary. For example, with two links crossing the boundary we can specify the holonomy between the corresponding two marked points on the boundary. This holonomy can be non-trivial even if the holonomy along the complete boundary is.
Remarkably, the result (5.7), applied to the BF vacuum with Z 2 gauge group, does agree with the entanglement entropy defined via the Hilbert space extension based on the spin network basis (or with the electric centre choice) found in [7]. To this end one has to choose m, the number of links connecting regions A and B to agree in both procedures. As we pointed out however, with our procedure, based on the fusion basis we are free to perform the (BF representation based) continuum limit keeping m fixed, ensuring a finite (or vanishing if m = 1) entanglement entropy for the BF vacuum. In contrast, using the Ashtekar-Lewandowski representation, the continuum limit of the BF vacuum state requires to take m to infinity, thus leading to an infinite entropy in this limit. Nevertheless this agreement in results is striking and it would be interesting to see if this holds for more generic states. In this case too, the contribution S A (D A (ρ ∂ )) vanishes as D A (ρ ∂ ) describes again a pure state. As shown in more detail in appendix C the other two contributions to the entanglement entropy are determined by the probability distribution Note that this agrees with the expectation value of the closed ribbon operator K[ρ ∂ ] along the boundary of the two regions, i.e. P (ρ ∂ ) = K[ρ ∂ ] . One has indeed ρ ∂ P (ρ ∂ ) = 1. With this probability distribution one can compute the entanglement entropy for the state under consideration to be More generally we can consider n non-intersecting ribbons R a [ρ a ], a = 1, . . . n, all going from region B to region A. These will generate a state proportional to the normalized state (see appendix C) We thus get the entropy to scale with the number of ribbon operators crossing the boundary. Note that we only get a non-vanishing entropy for these states due to the first two "classical" contributions to the entanglement entropy (5.1).
The distillable entanglement entropy would therefore be vanishing. This is probably due to the fact that the ribbon operators generating these states are not part of the reduced observable algebra O red , that underlies the definition of distillable entanglement entropy.

C. Comparison with the literature
A similar result to (5.11) was also obtained in [43,44], which considered the entanglement entropy for Chern-Simons theory. Whereas [43] uses the replica trick in a covariant path integral approach combined with surgery techniques, [44] employs again the replica trick but within the conformal field theory induced on the boundary.
In these two references, the entanglement entropy is also computed for states on the sphere as generated by the insertion of Wilson line operators. Notice, that the Chern-Simons Wilson lines involve holonomies of a Poisson noncommutative connection, which is therefore not the connection involved in our underlying states. At the same time, our states are generated by the action of ribbon operators, which are exactly Wilson lines for the double of the group. The analysis of [45,46] actually shows explicitly that the quantization of BF theory with gauge group G-like the one we are considering in this paper-is equivalent to the combinatorial quantization of Chern-Simons theory for the double D(G). This result 13 suggests that by appropriately identifying the structures in our computation and in theirs, they should match exactly, in the sense that they aim to compute the same physical quantities in two equivalent quantization schemes. After this premise, we can now compare the our results.
First of all, while [43] finds finite answers, those of [44]-found through computations in the edge field theory-are divergent. This divergence is due to an offset proportional to the central charge of the dual conformal field theory, and can therefore be thought as being associated to its zero-point energy. We will come back to this term later, for the moment let us focus on the other terms on which [43] and [44] agree. These terms contain two contributions. The first contribution S 1 coincides with our result (5.11), while the second S 2 is another universal offset determined by the total quantum dimension Ω S 2 = − ln Ω. (5.12) This can be interpreted as a 'vacuum contribution' to the entanglement entropy. It is a negative contribution. As such, it cannot result from a finite Hilbert space computation, and in fact we do not find it. (This term is also related to the so-called topological entropy [47,48], whose study in the framework of [48] needs the application of the fusion basis technique to disconnected regions. We postpone this study to future investigations.) Finally, coming back to the (positive) divergent term found using the dual field theory, it is interesting to speculate about its possible origin in the discrete framework. In Chern-Simons theory on a three-manifold with boundaries, the dual CFT living on the boundary is given by the WZNW model. If the three-manifold is a two-disk times an interval, a way to see the appearance of the WZNW model is by solving explicitly the flatness constraint in the bulk of the disk. In the standard notation A µ (x) = g(x) −1 ∂ µ g(x). This leaves us with the field g(x) on the boundary, since its bulk contributions to the action essentially cancel out (modulo topological terms). Therefore, it is the choice of local frame g(x) in which the flat gauge connection is evaluated which encodes the boundary CFT field (see also [49] for these derivations in relation to three-dimensional gravity). At this point, it is natural to draw a parallel between the frame g(x) in the continuum theory and the local frames we introduce in the refined boundary-puncture picture. If this is done, it is clear that infinitely many refining points on the boundary-puncture are needed to fully capture the dual field theory. In this limit also our entropy diverges. But as our procedure and the regularization procedure used in [44] are completely different, a more precise relation is difficult to obtain at this stage. In spite of this, we find this question extremely interesting.

D. TQFT based continuum limits
In section IV, we defined the Hilbert spaces H p that capture 2p − 2 degrees of freedom of a gauge theory on a fixed graph or lattice embedded on the 2-sphere. The set-up taken here allows to embed these Hilbert spaces H p into a continuum Hilbert space H cont .
This leads to the so-called BF representation [16][17][18], consisting of a Hilbert space which supports a representation of a continuum observable algebra, formed by the ribbon operators. The Hilbert space is based on the BF vacuum state, which is sharply peaked on flat connections, and is spanned by states that arise from the action of finitely many open ribbon operators on this vacuum state. The ribbon operators are then allowed to end at arbitrary points, which hence define the punctures. Such a Hilbert space can be constructed as an inductive limit from a family of Hilbert spaces based on fixed graphs (or more precisely equivalence classes of graphs, see [16][17][18]). In this latter viewpoint one puts all degrees of freedom, which are finer than the ones supported by the fixed graph, in the BF-vacuum state.
This leads us to the following interpretation of the result (5.11). First of all, the vacuum state has a vanishing entanglement entropy. Then, each (charge ribbon) operator R[ρ] that connects the two regions contributes to the entanglement entropy with ln dim ρ. This is despite the fact that in our definition of entanglement entropy we make explicit use of a particular fusion basis, which involves only one ribbon crossing the boundary.
An analogous result holds for the 'electric centre' choice, or its spin-network based extension [8]. Here we can also introduce a continuum Hilbert space, known as Ashtekar-Lewandowski representation [13][14][15]. This Hilbert space is based on a vacuum state peaked on vanishing electric fluxes. Wilson loops and lines do now act as creation operators. In fact a spin-network state results from the action of a network of Wilson lines connected via intertwiners. Using the spin-network based extension, one also finds that the entanglement entropy of a spin-network basis state is given by a ln dim ρ a where now ρ a denotes the irreducible representation of G associated to the a-th spin-network link which crosses the boundary. Again, one finds that the associated vacuum (the Ashtekar-Lewandowski vaucuum) has vanishing entanglement entropy.
However, if we express the BF vacuum in the spin-network basis and use the procedure of [8] to define the entanglement entropy, we notice that the result depends on the underlying graph. In particular, taking a refinement limit for the graph we would find a divergent result. Such a refining limit is necessary in this case to fully (i.e. everywhere) describe the BF vacuum. In other words, the BF vacuum is an infinitely excited state with respect to the Ashtekar-Lewandowski one, and the entanglement entropy reflects this fact. On the other hand, using the fusion basis and the related Hilbert space extension we emphasize the excitations relative to the BF vacuum itself. This leads to a result which is graph independent, a fact that makes our method applicable to the case of (2 + 1) dimensional gravity, which is described via a BF theory with defects describing point particles. Of course, attempting a description of the Ashtekar-Lewandowski vacuum (that is the strong coupling limit of Yang Mills theory) in terms of the fusion basis would also lead to results which are graph dependent or divergent.
Thus, we see that different notions of entanglement entropy are also adjusted to different notions of representations, or phases, or regimes. The BF representation corresponds to the (Yang Mills) weak coupling regime, and the Ashtekar-Lewandowski representation to the strong coupling regime. The excitations are in both cases (quasi-local) deviations from the weak coupling and strong coupling limit, respectively. In both cases, one deals with a topological theory with defect excitations.
We wish to emphasize that the vacua, both in the BF as well as in the Ashtekar-Lewandowski representation, describe theories without propagating degrees of freedom. That is, one can define Hamiltonians for which these vacua are the lowest energy states, which are moreover gapped. It is, thus, consistent to associate a vanishing entanglement entropy to these states. To describe the vacua of theories with propagating degrees of freedom in the continuum limit, we would need to introduce infinitely many excitations with respect to either of the vacua. This would lead to the usual divergent behaviour for the entanglement entropy in quantum field theories with propagating degrees of freedom.

VI. ENTANGLEMENT ENTROPY IN GRAVITY
The notion of entanglement entropy is usually (but not exclusively) associated to subsystems describing a region of space, which is specified by coordinates. One would like to link such a choice of region to a subset of observables, commuting with the remaining observables, but we have seen that this is already an ambiguous process for gauge systems. These difficulties are much more enhanced in general relativity.
In background independent theories, such as general relativity, regions specified by coordinates lack an a priori operational meaning. Alternatively, one can define regions through matter or metric fields. This is very similar to employing relational observables [50][51][52][53] as gauge invariant observables in general relativity. Here the metric or matter fields are used as a reference system, in which other fields can be expressed in [53][54][55]. Thus one can also attempt to specify 'physical regions' by employing a 'physical reference system'.
Relational observables can be computed in an approximation scheme [56,57], which also allows an understanding of how the standard observables of quantum field theory (on a fixed background) arise as approximations to fully gauge invariant observables.
A crucial drawback of using matter or metric fields as reference system is that there will be phase space regions in which these fields are not suited as clocks and rods. In some systems smooth gauge invariant observables might not exist [58,59]. Thus one expects that notions of locality can be realized only for a certain class of states [4] and are furthermore only approximate [29,51].
Another key point is the question whether one can find a split of the observable algebra into mutually commuting sets describing (approximately local) subsystems [4,29,60]. For example, the approximation scheme developed in [56] regains the usual quantum field theoretical observables on a fixed background at lowest order, but at higher orders it includes non-local terms. Giddings and Donnelly argue that observables creating e.g. matter fields, need to be gravitationally dressed, in order to capture the accompanying gravitational field. In contrast to Yang-Mills theories, this dressing cannot be screened and leads to an inherent non-local structure [29].
Using relational observables, one can deduce the commutator algebra by using Dirac brackets [52]. Realistic (that is relativistic) matter fields allow, however, only for an approximate localization [51,56]. Physical coordinates built from geometry (e.g. by using geodesics) lead, at least so far, to non-local algebras [56,61].
The exploration of the diffeomorphism invariant observable algebra is very difficult, as it basically requires to understand and solve the dynamics of the system. This is, of course, a very challenging task for the four-dimensional theory. On the other hand three-dimensional general relativity is much simpler: it describes locally flat spacetimes (or homogeneously curved ones, in presence of a non-vanishing cosmological constant). This also means that one has no local degrees of freedom, but only global topological ones. Introducing matter changes this situation, but it requires again a solution of the theory. To keep the system solvable, one can consider the coupling of point particles, which leads to a topological field theory with curvature and torsion defects, as discussed in this paper.
In fact (Euclidean) 3D gravity without a cosmological constant can be described by a BF theory with SU(2) structure group. The coupling of point particles leads to curvature and torsion defects [41,[62][63][64]. In short, the formalism needed to describe 3D gravity is very close to the formalism used here. Moreover including a positive cosmological constant, one has to work with a SU(2) q structure group, with q a root of unity. This leads to a finite dimensional Hilbert space, as in the finite group case we discussed. The fusion basis, ribbon operators, as well as gluing and cutting procedures are also available in the quantum group case [65].
The fusion basis diagonalizes a maximally commuting subset of gauge invariant (i.e. Dirac) observables, given by closed ribbon operators. A conjugated set of Dirac observables is provided by open ribbon operators going from one particle to another. Adopting the definition of entanglement entropy laid out here, a region is indeed specified by its matter content, that is by the particles contained in this region. Note that we do not have to specify the precise (geometric) position of the boundary, we only need to declare which particles belong to which regions. The geometric information is rather contained in the state under consideration.
We will refer to particles in 'region' A as A-particles and the remaining particles as B particles. The associated Aobservable algebra includes closed ribbons surrounding subsets of A-particles, with the exception of the closed ribbon surrounding all A-particles (and therefore also surrounding all the B particles as we consider spherical topology). This closed ribbon forms the centre of the algebra in the language of [9]. Additionally one can construct Dirac observables from open ribbon operators going from one A-particle to another A-particle.
Ribbon operators crossing the boundary cross also the closed ribbon along the boundary and would therefore not commute with it. These operators cannot be associated to either the A or B region.
Thus, although we can solve in this example the problem of how to define a region, we still have to modify the observable algebra, removing ribbons that cross from the A to the B regions from the operator algebra on which the entanglement entropy is being defined (if we follow the definition of [9]).
The notion of subsystems for 3D gravity used here differs in key points from the proposal (so far on the classical level) of [11]. There one introduces additional fields, that allow to fix the boundary in terms of embedding or coordinate functions. This has been motivated as a generalization of the extended Hilbert space construction (or rather its classical version). Here we point out that as there are different extension procedures in lattice gauge theories, this is also very likely to hold for gravity. The procedure laid out in this work can be applied to 3D gravity, and leads to much less extra structure compared to [11]. Furthermore, we can also state the definition of subsystems in terms of mutually commuting subsets of the Dirac observable algebra. This has still to be addressed within the proposal of [11], as was also remarked in [29]. Perhaps, the relevant suggestion is that this splitting of the observable algebra might not be only achieved by removing certain observables, but also allowing for (many) more observables via the introduction of a new unphysical 'boundary field'. This additional boundary structure might, in fact, be used to construct new local observables which otherwise would not be available [29].

VII. DISCUSSION
Recent work has shown that the notion of entanglement entropy in gauge systems is ambiguous. The deep underlying reason is that due to the non-local features of the observable algebra in gauge systems, a notion of subsystems needs to be defined first. The way this question is answered does not only affect the definition of entanglement entropy, but has much wider implications for our understanding of (quantum) systems with gauge symmetries [4,11,30,66]. In particular, defining subsystems in background independent theories, e.g. gravity, leads to various completely open issues. The methods developed here lead to a new proposal for lattice gauge theories, that is also applicable to (2 + 1) dimensional gravity. The main feature of this proposal is to use defect excitations to localize regions, which in the case of (2 + 1) gravity means that regions are specified operationally by their particle content.
Furthermore, we clarified the relation between the different approaches put forward so far to define entanglement entropy, notably the extended Hilbert space approach [8], and the CHR approach [9], which focuses on the observable algebra. In particular, we showed that the extended Hilbert space approach can be generalized to match not only the 'electric centre' choice of [9] but also the 'magnetic centre' choice [9] (and its non-Abelian generalization). In our view, the resulting notion for subsystems can, in both approaches, be fully characterized by a choice of boundary conditions. In the non-Abelian case, the extended Hilbert space approach relies on the introduction of extra frame information at the boundary, which we argued could also be added in a generalized CHR approach. We have also seen that the proposal made here requires only the introduction of a global frame, which is then transported with a locally flat connection along the boundary. In contrast, the spin-network based method of [8] necessarily introduces for each link cut by the boundary-and in the continuum limit to each point of the boundary-a local frame. Nevertheless, we observed that-if we wanted to-we could extend our framework as well, by allowing arbitrarily many frames along the boundary, leading to additional contributions to the entanglement entropy.
We have also pointed out that the different choices of boundary conditions, described by the 'electric' vs. 'magnetic' centre, are related to a choice of vacuum state. These vacuum states are of a 'topological' nature, i.e. they arise as vacua of topological field theories with no local degrees of freedom. The states can be used to define continuum Hilbert spaces, that then describe the states of the related topological field theory with defect excitations. Thus, for states describing BF theory with defects, we have to choose the (generalized) 'magnetic centre' definition, in order to obtain an entanglement entropy which is (i) regularization (i.e. graph) independent and (ii) finite.
The vacua we discussed here, are of a squeezed nature, which means they are sharply peaked either on flat connection (for the BF representation [16][17][18] ) or vanishing electric fluxes (for the AL representation [13][14][15]). The relation to preferred boundary conditions arises for the following simple reason: for states sharply peaked on connection degrees of freedom, it is natural and appropriate to fix the connection degrees of freedom (or the curvature) at the boundary. Similarly, for states peaked on some value of the electric flux, the original extended Hilbert space procedure [8] based on spin-networks is the most natural and appropriate one. Possible generalizations include q-deformed BF theory vacua [65], corresponding to (2 + 1) gravity with a cosmological constant and, in condensed matter, to string net models [48]. Furthermore, we suspect that also vacua with non-vanishing background values-e.g. for the electric fluxes [67][68][69]-come with a preferred notion of entanglement entropy.
Our methods employ techniques from topological field theories, which in their 'extended' form can be defined on manifolds with boundaries and corners. We believe that this direction can be further explored in order to learn how to define the notion of subsystems, in particular for background independent systems.
where it is understood that the symbols C ρ1ρ2ρ3 I 1 I 2 I 3 represent matrices C ρ1ρ2 whose indices are given by the composed labels ρ 3 I 3 and I 1 I 2 . By analogy with the group case, such maps are referred to as Clebsch-Gordan coefficients. From the unitarity of C ρ1ρ2 , it follows the orthogonality relation I1,I2 C ρ1ρ2ρ I1I2I · C ρ1ρ2ρ I1I2I = δ ρ,ρ δ I,I , as well as the completeness relation ρ I C ρ1ρ2ρ I 1 I 2 I · C ρ1ρ2ρ I1I2I = δ I 1 ,I1 δ I 2 ,I2 .