Geometric Entanglement in Topologically Ordered States

Here we investigate the connection between topological order and the geometric entanglement, as measured by the logarithm of the overlap between a given state and its closest product state of blocks. We do this for a variety of topologically-ordered systems such as the toric code, double semion, color code, and quantum double models. As happens for the entanglement entropy, we find that for sufficiently large block sizes the geometric entanglement is, up to possible sub-leading corrections, the sum of two contributions: a bulk contribution obeying a boundary law times the number of blocks, and a contribution quantifying the underlying pattern of long-range entanglement of the topologically-ordered state. This topological contribution is also present in the case of single-spin blocks in most cases, and constitutes an alternative characterisation of topological order for these quantum states based on a multipartite entanglement measure. In particular, we see that the topological term for the 2D color code is twice as much as the one for the toric code, in accordance with recent renormalization group arguments. Motivated by this, we also derive a general formalism to obtain upper- and lower-bounds to the geometric entanglement of states with a non-Abelian symmetry, and which we use to analyse quantum double models. Furthermore, we also analyse the robustness of the topological contribution using renormalization and perturbation theory arguments, as well as a numerical estimation for small systems. Some of our results rely on the ability to disentangle single sites from the system, which is always possible in our framework. Aditionally, we relate our results to the relative entropy of entanglement in topological systems, and discuss some tensor network numerical approaches that could be used to extract the topological contribution for large systems beyond exactly-solvable models.

Here we investigate the connection between topological order and the geometric entanglement, as measured by the logarithm of the overlap between a given state and its closest product state of blocks. We do this for a variety of topologically-ordered systems such as the toric code, double semion, color code, and quantum double models. As happens for the entanglement entropy, we find that for sufficiently large block sizes the geometric entanglement is, up to possible sub-leading corrections, the sum of two contributions: a bulk contribution obeying a boundary law times the number of blocks, and a contribution quantifying the underlying pattern of long-range entanglement of the topologically-ordered state. This topological contribution is also present in the case of singlespin blocks in most cases, and constitutes an alternative characterisation of topological order for these quantum states based on a multipartite entanglement measure. In particular, we see that the topological term for the 2D color code is twice as much as the one for the toric code, in accordance with recent renormalization group arguments [H. Bombin, G. Duclos-Cianci, D. Poulin, New J. Phys. 14 (2012) 073048]. Motivated by these results, we also derive a general formalism to obtain upperand lower-bounds to the geometric entanglement of states with a non-Abelian group symmetry, and which we explicitly use to analyse quantum double models. Furthermore, we also provide an analysis of the robustness of the topological contribution in terms of renormalization and perturbation theory arguments, as well as a numerical estimation for small systems. Some of the results in this paper rely on the ability to disentangle single sites from the quantum state, which is always possible for the systems that we consider. Additionally we relate our results to the behaviour of the relative entropy of entanglement in topologically-ordered systems, and discuss a number of numerical approaches based on tensor networks that could be employed to extract this topological contribution for large systems beyond exactly-solvable models.

I. INTRODUCTION
Topological order (TO) [1] is an example of new physics beyond Landau's symmetry-breaking paradigm of phase transitions. Systems exhibiting this new kind of order are linked to concepts of the deepest physical interest, e.g. quasiparticle anyonic statistics and topological quantum computation [2]. Importantly, TO finds a realisation in terms of topological quantum field theories [3], which are the low-energy limit of quantum lattice models such as the toric code and quantum double models [4], as well as string-net models [5].
A remarkable property about TO is that it influences the long-range entanglement in the wave function of the system. For instance, as proven for systems in two spatial dimensions (2D) [6,7], the entanglement entropy S of a region, such as A in Fig. 1.a, of boundary size L 1 obeys the law where S 0 ∝ L is a non-universal term (the so-called "boundary law", or "area law"), ν some exponent, and S γ is a universal long-distance contribution: the topological entanglement entropy. S γ is non-zero for systems with TO, e.g. S γ = 1 for systems in the topological phase of the toric code. In general one has that where {d a } are the so-called quantum dimensions of the associated anyon model. More recently, similar universal contributions have also been found for other bipartite entanglement measures such as the mutual information and the Rényi entropy [8,9]. The main purpose of this paper is to investigate, in the context of systems with TO, a global measure of entanglement which captures multipartite correlations in the system: the geometric entanglement (GE) [10,11], which we shall denote by E G . This measure intuitively characterises how well an entangled state can be approximated by a mean-field state, i.e., a product state, and its behaviour is in principle very different from that of bipartite measures such as the entanglement entropy. Can such a mean-field, global measure of entanglement be able to reveal any property about topological order? As we shall see in this paper, this is indeed the case.
Most of our study is focused on exactly-solvable models corresponding to fixed-points of the renormalization group (RG). We deal most prominently with the toric code model [4], but shall also discuss the double semion model [5], topological color codes [12], as well as quantum double models [4]. Each one of these models can be recast as some particular instance of a string-net (Levin-Wen) model. Because of this, they correspond (by definition) to RG fixed-points and, therefore, they are representative of their respective topological phases. Our main result is that, for all fixed-point models considered in this paper, the GE obeys with E γ being identical to S γ and E 0 some bulk contribution. The exact form of E 0 depends on short-distance details such as the size of the blocks (i.e. the number of spins) considered for the closest product state optimisation. As we shall see, for blocks of boundary size L, E 0 will be proportional to L and also to the number of blocks n b (more specifically, E 0 ∝ n b L). We call this behaviour, which depends on short-distance correlations, a "boundary law" for the GE, since it corresponds to a boundary term for each one of the blocks, i.e., E 0 /n b ∼ L.
Moreover, E γ being identical to S γ is an indication of the topological origin of the term E γ . Although we do not have a proof that Eq. (3) holds generically beyond the considered fixed-point states, we shall argue that the above law is also valid up to subleading corrections away from the RG fixed points, e.g. for the toric code model under the influence of nonrelevant (and weak) perturbations such as magnetic fields. In this case the geometric entanglement of blocks with boundary size L 1 seems to obey the law where again E 0 ∝ n b L is a bulk contribution, ν some exponent, and E γ the topological contribution. From here on, we shall refer to the term E γ as the topological geometric entanglement. It is worth mentioning now a number of remarks about our results. First, they provide the first evidence that topological order can actually be read and assessed in a number of situations from multipartite nature of entanglement, and, in particular, the geometric measure of entanglement employed here. Thus, this may provide an alternative way of characterising topological order in quantum states. Second, our results are also the first explicit and analytic examples of a boundary law behaviour for the GE in ground states of 2D quantum many-body systems. Third, we will see that the topological term for topological color codes is twice the one for the toric code. This is consistent with the recent result that the former model is equivalent to two copies of the latter [13]. Moreover, there may be potential advantages in considering the GE to characterise topological order instead of other quantities. For instance, one of them is that given a (possibly approximate) description of the ground state by a Tensor Network such as a Projected Entangled Pair State (PEPS) [14] or a Multi-scale Entanglement Renormalization Ansatz (MERA) [15], its computation is not too costly. This is because, in the context of Tensor Network methods, the calculation mainly involves overlaps with product states and optimisations over these, which can be computed quite efficiently. In this respect, some of the possible approaches using these methods will be discussed in the Appendix. Moreover, the GE can also be probed experimentally with current technology using e.g. Nuclear Magnetic Resonance via single spin measurement on all nuclear spins [16], or ultra-cold atoms in optical lattices via similar techniques to those explained in Ref. [17].
The structure of the paper is as follows. First, we briefly review the basics of the geometric entanglement in Sec. II. After this, we deal in Sec. III with the toric code model. Since this is the simplest model displaying non-trivial topological order, we spend quite some time explaining many of its properties, as well as derivations of the geometric entanglement of spins and blocks both for the square and honeycomb lattices. In Sec. IV we explain a number of important points that need to be considered in order to generalise the results obtained for the toric code to other models. Then, the double semion model is analysed in Sec. V, and topological color codes are addressed in Sec. VI. For color codes we derive results using two independent methods: a bound approach (similar to the one used for the other models), and a direct approach valid for Calderbank-Shor-Steane (CSS) self-orthogonal codes. Depending on the particular setting we will see that one approach may provide some advantages over the other. After this, we develop in Sec. VII a general formalism for the bound approach, valid for states with non-Abelian symmetries, and apply this formalism to quantum double models in Sec. VIII. The situation beyond RG fixed points is considered in Sec. IX, where the robustness under perturbations of the topological contribution is assessed by RG and perturbation theory arguments, as well as by small-size numerical calculations. We present our conclusions and possible future directions in Sec. X. Finally, in the Appendix we describe how our results are related to those of the relative entropy of entanglement for topologically-ordered states, and discuss possible numerical algorithms based on tensor networks to extract this topological component for non-exactly solvable models of large size.

II. GEOMETRIC ENTANGLEMENT
To begin with, let us recall some basic notions about the geometric entanglement. Consider an m-partite nor- is the Hilbert space of party i. For instance, in a system of n spins each party could be a single spin, so that m = n. But each party could also be a set of spins, either contiguous (a block [18]) or not.
We wish now to determine how well state |Ψ can be approximated by an unentangled (normalized) state of the parties, |Φ ≡ ⊗ m i=1 |φ [i] . The proximity of |Ψ to |Φ is captured by their overlap. The entanglement of Bipartition with region A and region B. Region A has linear size L. This is usually the partition used in evaluations of the entanglement entropy and other bipartite entanglement measures. (b) Multipartition with many blocks, each one with linear size L. This is the partition used in the calculations of the geometric entanglement. If each block contains only one spin (akin to individual degree of freedom), then the product state is a completely separable state and the geometric entanglement is defined with respect to individual spins. However, each block can also contain multiple spins, and in this case the geometric entanglement is defined with respect to the closest product state of the individual blocks. In this latter scenario, we can inquire how the entanglement depends on the block boundary size L and, in particular, whether there is something like a boundary law (i.e. linear scaling with L) and whether there is a topological correction.
|Ψ is thus revealed by the maximal overlap [10], The larger Λ max is, the less entangled is |Ψ . Therefore, we quantify the entanglement of |Ψ via the quantity where we have taken the base-2 logarithm, and which gives zero for unentangled states. E G (Ψ) is called geometric entanglement (GE). This quantity has been studied in a variety of contexts, including critical systems and quantum phase transitions [18,19], quantification of entanglement as a resource for quantum computation [20], local state discrimination [21], and has been recently measured in NMR experiments [16].
To make a connection between the GE and other measures of entanglement, let us consider the case of just two sets of spins. In this case E G (Ψ) coincides with the so-called single-copy entanglement between the two sets, with ν 1 (ρ) the largest eigenvalue of the reduced density matrix ρ of either set [22]. As is well known, this also coincides with the α-Rényi entropy for α → ∞.
Unlike other measures of entanglement, the GE offers a lot of flexibility to study multipartite quantum correlations in spin systems. For instance, one can choose each party to be a single spin, but one can also choose blocks of increasing boundary length L [18] (see Fig. 1.b for an example). In fact, studying how the GE changes with L provides valuable information about how close the system is to a product state under coarse-graining RG transformations.

III. TORIC CODE MODEL
As described in the introduction, here we considered a variety of models corresponding to different topological phases. The simplest of these models is the toric code [4]. This model is, in fact, the simplest example of a topologically non-trivial 2D system, and our results for this model provide also the grounds for the rest of the models that we shall consider (namely double semion, color codes, and quantum double models). Given its relevance for the rest of the paper, we provide now a short introduction to some of the key aspects of this model. As we shall see, this will be very useful for the forthcoming calculations.
Mathematically speaking, the toric code is the RG fixed point of the topological phase of a Z 2 gauge theory. The model is equivalent under local transformations (or local moves, or local disentanglers) to the Levin-Wen string model on a honeycomb lattice [5], which is by definition a RG fixed point. A more precise derivation of this property will be provided later.
To define the toric code, we consider a lattice Σ on a torus. In this paper we will deal with the toric code in the square and honeycomb lattices, yet the model can be defined on arbitrary lattices. Other Riemann surfaces of genus g could also be considered easily without changing our conclusions. There are spin-1/2 (qubits) degrees of freedom attached to each link in lattice Σ. The model is described in terms of stars and plaquettes. A star "s" is a set of links sharing a common vertex. A plaquette "p" is an elementary face on the lattice Σ. For any star s and plaquette p, we consider the star operators A s and plaquette operators B p defined as where σ [j] α is the α-th Pauli matrix at link j of the lattice. Let us call respectively n s , n p and n the number of stars (vertices), plaquettes (faces) and links (sites) in lattice Σ. Importantly, star and plaquette operators satisfy the global constraint Therefore, there are n s − 1 independent star operators and n p − 1 independent plaquette operators.
With the definitions above, the Hamiltonian of the model reads This Hamiltonian is frustration free (i.e. all the terms in the above sum commute with each other), and can be diagonalized exactly as explained in Ref. [4]. The ground level is 4-fold degenerate on a torus (4g-fold for a Riemann surface of genus g). This degeneracy depends on the underlying topology of lattice Σ, and is already by itself a signature of topological order. Moreover, the ground level is a stabilised space of G s , the group of all the possible products of independent star operators, of size |G s | = 2 (ns−1) .

A. Ground states
In order to build a basis for the ground level subspace, let us consider a closed curve γ running on the links of the dual lattice Σ * . We define its associated loop operator x , where j ∈ γ are the links in Σ crossed by the curve γ connecting the centres of the plaquettes. It is not difficult to see that the group G s can also be understood as the group generated by all the possible contractible loop operators [41]. Let also γ 1 and γ 2 be the two non-contractible loops on a torus, and define the associated string operators w 1,2 ≡ W x [γ 1,2 ]. We call |0 and |1 the eigenstates of σ z respectively with +1 and −1 eigenvalue. With this assumptions, the ground level subspace then reads L = span{|i, j , i, j = 0, 1}, where It is easy to check that the four vectors |i, j are orthonormal and stabilised by G s , that is, g|i, j = |i, j ∀g ∈ G s and ∀i, j. These four states form a possible basis of the ground level subspace of the toric code model on a torus.

B. Excited states
Excited states of the toric code can be constructed by locally applying Pauli operators σ x , σ y , σ z on the ground states |i, j . Pauli operators σ z create pairs of deconfined charge-anticharge quasiparticle excitations, σ x create pairs of deconfined flux-antiflux quasiparticles, and σ y creates a flux-antiflux and a charge-anticharge pairs. These operators can be applied to several sites of the lattice, thus creating a quasiparticle pattern that defines the excited state. The excited states of the model are then labeled by a set of quantum numbers, |φ, c, i, j , where φ and c are patterns denoting the position of flux-type and charge-type excitations respectively, and i, j label the ground state |i, j that was excited. In this paper we will focus on the entanglement properties of states |φ, c, i, j . CNOT operations that disentangle qubits from the system, thus removing their links from the lattice. In the diagram, the arrows go from controlling to target qubits; (c) Example of a partition into blocks of a given boundary L, e.g. L = 12 here (but all our derivations work for arbitrary L). Qubits crossed by the boundary of a block are regarded as being inside the block; (d,e) Doing CNOT operations locally inside of each block we can remove all the stars inside of all the blocks in two steps: first, we apply CNOTs as in (a) to get (d) and, second, we apply CNOTs as in (b) to get (e). Notice that in (e) each star is made of either 4 or 8 qubits. The procedure works similarly for other blockings.

C. Disentangling the toric code
A key property of the toric code model, which is fundamental for some of the derivations in this paper, is that its ground state |0, 0 ≡ |i = 0, j = 0 can be created by a quantum circuit that applies a sequence of Controlled-NOT (CNOT) unitary operations over an initial separable state of all the qubits [23]. This means that it is actually possible to disentangle qubits from the ground state of the system simply by reversing the action of these CNOTs. Moreover, these CNOT operations are local, i.e. they do not span over delocalized sites in the lattice, but rather act over pairs of not-too-far-neighbouring qubits (e.g. nearest-and next-to-nearest-neighbours). This was the key observation that allowed to build an exact MERA representation of the ground states |i, j of the model [24]. The two fundamental disentangling moves are represented in Fig 2.a-b, and leave the overall quantum state as a product state of the disentangled qubits with the rest of the system. What is more, the rest of the system is left in the ground state of a toric code model on a deformed lattice Σ, where Σ is obtained from Σ by removing the links that correspond to the disentangled qubits. This property turns out to be of great importance for some of our derivations. , and back to a square lattice Σ . Each red arrow represents a CNOT operation that disentangles qubits from the system, as explained in Fig. 2. The arrows go from controlling to target qubits.

D. GE of the toric code
We are now in position to study the GE of the toric code model. Given that this is a paradigmatic model of topological order, we will perform a detailed analysis. First, we will provide a number of lemmas and theorems that will bring us towards an expression for the GE of the toric code in the square lattice, both for spins and blocks, in a particular basis. Second, we will make a similar study to the case of the honeycomb lattice. There we will see that different choices of ground states may give rise to different expressions for the GE of spins with the block size being unity, yet the long-wavelength properties are the same and equal to those for the square lattice after blocking of spins are considered. This is expected, as both lattices can be mapped between each other by means of CNOT disentangling operations (as shown in Fig. 3), and therefore they should correspond to the same RG fixed point. Here, we would like to emphasize that blocking is essential as the topological contribution is a long-wavelength property, and, thus, that the behaviour of the GE for sufficiently large block sizes is the same independently of which of the two lattices we choose. (Without blocking, this conclusion cannot be reached.)

A couple of Lemmas
Let us start by considering two Lemmas that will be quite useful in our derivations:  Proof. Let us choose a reference state from the ground level basis, e.g. |0, 0 . First of all, notice that the four basis states |i, j in the ground level are related to |0, 0 by local unitary operators, since the string operators w 1 and w 2 are tensor products of 2 × 2 identity operators and Pauli-x matrices. Since these operators act locally on each spin, they do not change the entanglement of the quantum state. Moreover, the excited states |φ, c, i, j are created locally by applying tensor products of 2 × 2 identity and Pauli-x and z operators to |i, j . The states resulting from these operations have then the entanglement properties of state |i, j , which in turn are the same as those of the reference state |0, 0 . Thus, all states |φ, c, i, j have the same entanglement properties.

Lemma 2.
Consider an arbitrary set of blocks of qubits in lattice Σ, where each block is regarded as an individual party. Then, the entanglement properties of the ground state |0, 0 are the same as those of state | 0, 0 , the ground state of the toric code model in the deformed lattice Σ obtained after disentangling as many qubits as possible using CNOTs locally inside of each block.
Proof. Given the ground state |0, 0 in Σ, we start by doing CNOT disentangling operations locally inside of each block in order to disentangle as many qubits as possible. As a result, state |0, 0 transforms into state |0, 0 disentangled ≡ |e 1 ⊗ · · · ⊗ |e p ⊗ | 0, 0 , where |e k is the quantum state for the k-th disentangled qubit, k = 1, . . . , p (p is the number of disentangled qubits), and | 0, 0 is the i = 0, j = 0 ground state of a toric code Hamiltonian in the deformed lattice Σ. Since all CNOT operations are done locally inside of each block, states |0, 0 and |0, 0 disentangled have the same entanglement content if we regard the blocks as individual parties. What is more, the entanglement in the system is entirely due to the qubits in | 0, 0 , which proves the lemma.

Square lattice
We are now in position to consider the entanglement properties of the four ground level states |i, j as well as the excited states |φ, c, i, j . Our main result here is that, for these states, the GE of blocks consists of a boundary term plus a topological contribution. Keeping the above two lemmas in mind, we now present the following theorem about the geometric entanglement of spins in the square-lattice toric code: Theorem 1. For the toric code Hamiltonian in a square lattice Σ, the four ground states |i, j for i, j = 0, 1 and also the excited states |φ, c, i, j have all the same GE of spins and is given by where n s is the number of stars in Σ.
Proof. Using Lemma 1 we can restrict our attention to the GE for the state |0, 0 . In this setting, we call computational basis the basis of the many-body Hilbert space constructed from the tensor products of the {|0 , |1 } local basis for every spin, which is an example of product basis for the spins. Our proof follows from upper and lower bounding the quantity E G (0, 0) in Eq.(6) for the ground state |0, 0 . First, from the expression in Eq. (12) for the ground states we immediately have that the absolute value of the overlap with any state of the computational basis is |G s | −1/2 , and therefore Λ max ≥ |G s | −1/2 . From here, we get E G (0, 0) ≤ log 2 |G s |, which gives an upper bound.
Next, to derive a lower bound we use the fact that if we group the n spins into two sets A and B, then means that a partition of the system with respect to the two sets A and B is considered. Thus, we have that The trick to find a useful lower bound is to find an appropriate choice of sets A and B. In our case, we consider e.g. the bipartition shown in Fig. 4.a for even × even lattices (other cases can be considered similarly). Then, we use a Lemma by Hamma, Ionicioiu and Zanardi in Ref. [6], that the reduced density matrix is the subgroup of G s acting trivially on subsystem A/B. Importantly, for A and B chosen as explained above, it happens that G s (A) and G s (B) are trivial groups consisting of only the identity element. With this in mind, we see that the reduced density matrix ρ A has eigenvalues either zero or |G s | −1 , which is |G s |-fold degenerate. This means that (Λ Combining the two bounds, we get E G (0, 0) = log 2 |G s | = n s − 1, and from here Eq.(13) for the GE for spins follows immediately.
Importantly, the techniques used in Theorem 1 can also be used to deal with blocks of spins whenever the blocks form a bipartite lattice. In general, one may first disentangle as many qubits as possible inside the blocks. Then the remaining qubits are left in an entangled state where the GE is equal to the number of remaining independent star operators, which amounts to a boundary law term for each block (possibly with a sub-leading correction depending on how the blocks are chosen) plus a topological term. As an example of this, let us present the following theorem: Theorem 2. Given a partition of the square lattice Σ into n b blocks of boundary size L as indicated in Fig. 2.c (where L is measured in number of qubits), then the four ground states |i, j for i, j = 0, 1 and also the excited states |φ, c, i, j of the toric code Hamiltonian have all the same GE of blocks and is given by Proof. First, notice that Lemma 1 and Lemma 2 imply that we can entirely focus on the ground state | 0, 0 of a toric code model in the deformed lattice Σ from Fig. 2.e. As shown in Fig. 2.d-e, it is always possible to remove all the stars inside of each block in Σ just by doing CNOTs locally inside of each block. As a result, the stars in the deformed lattice Σ correspond to those in Σ that lay among the blocks. It is easy to see that, for n b blocks of boundary L, there are n s = n b L/4 of such stars. The rest of the proof follows as in Theorem 1, by upper and lower bounding E G (0, 0) for the partition with respect to blocks. For the upper bound, we use the fact that Λ [blocks] max ≥ Λ [spins] max ≥ | G s | −1/2 , where the first inequality again used the property that if larger blocks are considered then the maximum overlap is also larger, and G s is the corresponding group of contractible loop operators on Σ. From here E G (0, 0) ≤ log 2 | G s | follows. For the lower bound, we divide the system into two sets A and B of blocks as indicated in Fig.4.b. Following a similar reasoning as in Theorem 1, it is possible to see again that no element g ∈ G s will act trivially on A or B except for the identity element. From this point, the rest of the proof is simply equivalent to the proof for Theorem 1.

Interpretation of the previous Theorems
Let us now discuss Eqs. (13) and (14) in detail. First, let us remind that a variety of works have shown the existence of a "boundary law" for the entanglement entropy of many 2D systems [25], including models with TO [6,7]. It is then remarkable that, according to the first term of these equations, the entanglement per block obeys also a boundary law, i.e: it is proportional to the size L of the boundary of the block. To the best of our knowledge, these results are the first example of a boundary law behaviour for a multipartite (rather than bipartite) measure of entanglement in 2D.
However, the second term in Eqs. (13) and (14) is far more intriguing and important. Its existence is caused by the global constraint from Eq. (10) on star operators which, in turn, allow for the topological degeneracy of the ground state. Thus, this term is of topological nature, and quantifies the pattern of long-range entanglement present in topologically ordered states.
In order to clarify further the meaning of the topological term, let us consider the bipartite case of a block of spins of boundary size L and its environment (i.e. the rest of spins). In the bipartite case, the geometric entanglement E G coincides with the single-copy entanglement E 1 = − log ν 1 (ρ), with ν 1 (ρ) the largest eigenvalue of the reduced density matrix ρ of the block. Since this density matrix is proportional to a projector (see Ref. [6] and also Ref. [9]), we have that E 1 = S, with S the entanglement entropy of the bipartition. Thus, for this case we have the chain of equalities We also know that for a system with TO, the entangle- , where S 0 is some boundary law term and S γ is the topological entropy. As explained in Ref. [6], for the toric code model the global constraints in Eq.(10) imply that S γ = 1. Thus, with the convention from Eq.(4), we have that the topological geometric entanglement is given, in this case, by The topological origin of E γ is thus clear. Notoriously, this topological contribution is maintained when promoting the entanglement measure from the bipartite to the multipartite scenario.

Honeycomb lattice
For the honeycomb lattice, most of the derivations are equivalent to the square lattice, yet there is a small difference: while the square lattice is self-dual, the honeycomb lattice is not. This means, among other things, that the toric code in the honeycomb lattice is not self-dual with respect to a change of star-and plaquette-operators (which is actually the case for the square lattice). In practice, this implies that different choices of ground state basis may have different values of the GE of spins. However, this is only a short-distance property, since the GE of blocks has the same topological contribution as for the square lattice. In order to illustrate this more explicitly, we provide calculations for two specific ground states.
The |0, 0 ground state.-We first consider the |0, 0 ground state obtained as which is nothing but the i = 0, j = 0 ground state in the notation of Eq. (12). This time, G s is the group generated by the product of all possible star operators in the honeycomb lattice. As before, |G s | = 2 ns−1 , where n s is the number of stars (vertices) in the honeycomb lattice.
For this state, we provide the following observation (not a Theorem): Observation 1. For the toric code Hamiltonian in a honeycomb lattice Σ, the GE of spins for the ground state |0, 0 is upper-and lower-bounded by where n p is the number of plaquettes (faces) in Σ. Moreover, numerically it looks like the upper bound is saturated, so that The above is an observation, and not a theorem, since Eq. (19) has only been checked numerically. Yet, let us prove now the upper and lower bounds in Eq. (18). First, the maximal overlap can be easily upper bounded again by Λ max ≥ |G s | −1/2 = 2 (1−ns)/2 . Nevertheless, unlike in the case of the square lattice, this time we can actually have a better bound, i.e., which is also smaller than |G s | −1/2 = 2 (1−ns)/2 using the fact that n = 3n p and n s = 2n p for the honeycomb lattice. This gives the desired upper bound on the GE, To obtain the lower bound on GE (or equivalently upper bound on the maximal overlap), our strategy is again similar to the one for the square lattice. This time we  consider the bipartition illustrated in Fig. 5. There is a star operator in each of the smallest units of the sets A (red) and B (blue). It is easy to convince oneself that the nontrivial subgroup G s (A) of G s is generated by the n s /4 star operators residing in all units in A, and similarly for G s (B). The size of these groups is therefore which leads to n p − 1 ≤ E G , as claimed. Combining this bound with the one computed previously, we arrive at the result that we mentioned before, namely, n p − 1 ≤ E G ≤ n p + 1.
In order to see whether any of these bounds is tight, we have carried out a numerical computation for small system sizes and confirmed that it is the upper bound which is saturated, i.e., E G = n p + 1. This can be seen in Table I. We would like to emphasize that even though we do not obtain −1, this is only a short-range property (where each block size is unity), as E G for blocks will give rise to the correct contribution; see below.
The |+, + ground state.-We now consider the |+, + ground state. By definition, this ground state is obtained by acting on the reference state |+ ⊗n with the elements of the group G p of plaquette operators, namely, Notice that, unlike in the square lattice (where there is a duality between stars and plaquettes), in the honeycomb lattice the |+, + ground state does not need to have, necessarily, the same properties as the |0, 0 ground state from the previous section (at least not the same shortdistance properties). For this state we can actually find its GE exactly, and provide the following Theorem:

|+, + is given by
where n p is the number of plaquettes (faces) in Σ.
Proof. We can obtain the GE analytically as follows: first, we use the duality that the sate |+, + in the honeycomb lattice Σ is equivalent to the |0, 0 ground state of the toric code on the dual lattice Σ * , which is the triangular lattice. After this, we simply follow the same approach as before by finding appropriate upper-and lowerbounds. The trivial upper bound is given by n p −1 ≤ E G , with n p the number of plaquettes in the original honeycomb lattice. To find the lower bound, we consider the bipartition given in Fig. 6, which is similar to the one that we used square lattice (Fig. 2). Using the same methods as before, the expression for the GE follows.
Blocking the honeycomb lattice.-The fact that the |0, 0 state on different lattices gives different constant contributions to the GE of spins turns out to be a shortdistance property that disappears when blocks of spins are considered. Blocking corresponds to coarse-graining and gives long-wavelength properties. This is easy to notice, if one realises that in just one RG step we can map the honeycomb to the square lattice. Let us explain this in more detail. It is well known that the toric code in the honeycomb and square lattices can be mapped to one another by simple RG (disentangling) operations acting locally on the lattice. These transformations (from the square to the honeycomb lattice, and back to the square) are represented in the diagram of Fig. 3. In fact, these moves can also be understood in terms of Entanglement Renormalization steps [24], where some of the initial qubits become disentangled from the rest after every step. Importantly, after mapping the honeycomb to the square lattice, the GE follows exactly the laws explained in the previous section, and hence it is evident that E γ = 1 also holds in this case after a suitable blocking of the spins (so that the RG CNOTs fall within the blocks).
For instance, one could choose the blocking shown in Fig. 7 and obtain E γ = 1, but similar results hold for other choices of blocks as well. One could also consider larger blocking such as the one shown in Fig. 8. It is interesting to know that one-step of blocking is sufficient to remove the short-range behaviour and result in a correct topological contribution. Needless to say, in the process of blocking the corresponding "boundary laws" for the bulk contribution of the GE of blocks also follow.

Brief discussion of the honeycomb lattice results
Let us comment briefly on the above results. First, we have seen that the GE of spins is, in this case, different for the |0, 0 and |+, + states. One could also infer from here that the topological contribution to the GE is also different. However, we wish to stress that the topological GE is a long-range contribution, and can only be extracted reliably after appropriate blocking of the spins. This is exactly what has been observed when considering the blocking of spins in the honeycomb lattice. Second, it is clear that the same conclusions that applied to excited states and other ground states on the square lattice apply also to the honeycomb lattice, i.e. once chosen a ground state basis, all the states in the basis as well as the corresponding (quasiparticle) excitations have the same amount of entanglement as the "reference" ground state (which in this case is either |0, 0 or |+, + ), as long as they are related by local operations.

A different topology: the sphere
In addition to the topology of a torus, we also consider that of a sphere. The simplest geometry is a tetrahedron, shown in Fig. 9a. The toric code ground state on this topology is unique.
Using the same technique of upper and lower bounding the maximal overlap, one can easily show that its GE of spins (i.e. one-site blocks) is given by E G = n s − 1, where n s = 4 is the number of vertices and is related to the number of spins n via n = d · n s /2, where d = 3 is the degree of a vertex. The bipartition that gives rise to this result is shown in Fig. 9b.
In order to obtain a scaling we add one tetrahedron to each face, which can be iterated in order to add more sites to the system. The first step of such a procedure is shown in Fig. 9c. In general it is easy to obtain E G = n s − 1 = 2n/d − 1, where n s , n and d are, respectively, the number of vertices, the number of spins (edges), and the vertex degree of the associated geometry. This result shows that the -1 contribution is not only present in the torus topology, but in the sphere topology as well.

IV. BEYOND THE TORIC CODE
After our detailed study on the toric code model, we now turn to study the GE in other topological models that also correspond to fixed points of RG. From now on, some parts of our study will not be as detailed as for the toric code, and will focus mainly on some of the fundamental properties. Yet, in order to make everything as clear as possible, we need to see how some of the main properties of the toric code apply also to other models. In particular, there are three main points to consider: 1. Excitations and other ground states.-As for the toric code, excited states and other ground states that can be obtained by local operations acting on some reference ground state will share the same entanglement properties (c.f. Lemma 1). Therefore, in many cases we shall only discuss the entanglement of some "reference" ground state.
2. Disentangling general stabiliser states.-Very importantly, the disentangling property is not exclusive of the toric code only, but it holds for any stabiliser state as well. More concretely, any state that is stabilised by a set of operators can be built from a quantum circuit acting on some separable initial state, where the gates of the quantum circuit involve the basis of the stabiliser group (see e.g. [26]). If this basis of stabilisers can be thought of as acting locally on a lattice, then the quantum gates of the corresponding quantum circuit will also be local in the lattice. Moreover, there will be a set of unitary gates in the quantum circuit that will correspond to each one of the elements of the group basis, and which will take into account the existence of that particular basis element. That is: removing the action of the gates corresponding to a particular basis element effectively removes that stabiliser from the basis of the stabiliser group (e.g. removing a plaquette or star operator from the toric code). Therefore, simply by reversing the action of such a quantum circuit on the stabilised state, we will always manage to disentangle the state by means of sufficiently local unitary operations. At every step, the remaining quantum state will be the state stabilised by the remaining set of stabiliser operators. In particular, this implies that there is always a MERA representation for any stabiliser state (not necessarily topological), and the toric code model is just a particular case.
This property turns out to be quite important, since it allows us to extend straightforwardly the results for the GE of blocks beyond the toric code: when blocking is considered, simply perform disentangling operations inside of each block in order to remove the maximum of stabilisers within each block, so that the remaining stabilisers will always be between spins in different blocks. By construction, the remaining state will usually obey an area law for the GE of blocks, plus a possible topological term. Because of this, from now on the GE of blocks will not be discussed in detail except for a couple of examples.
3. The non-Abelian case.-It is indeed possible to generalise our results to systems with non-Abelian symmetries, such as quantum double models. However, in order to deal with these models, it will be worthwhile to develop part of our formalism in the most general fashion. This will be done in Sec. VII, and quantum double models will be considered subsequently in Sec. VIII.

V. DOUBLE SEMION MODEL
The double semion model is given by the spin model on the honeycomb lattice where A s and B p are mutually commuting and given by (25) This is, the star operators are the same as for the toric code, whereas the plaquette operators are related via with B p the toric code plaquette operator. As for the toric code, star and plaquette operators satisfy the nonlocal constraint This model is known to be also topologically ordered, yet its topological phase differs from the one of the toric code model. More precisely, while the toric code corresponds to the topological phase of a Z 2 gauge theory, the double semion model corresponds to a U (1)×U (1) Chern-Simons theory [5].
One of the ground states of this double semion model is almost identical to the toric code ground state |+, + , except that each term in the decomposition Eq. (22) is multiplied by a factor (−1) Xc , where X c counts the number of loops formed by spins in the |− state. More specifically, this ground state is given by where again G p is the group of all possible products of plaquette operators B p . Of course, the state |+, + is by construction a stabilised state of the group G p . We provide now the following Theorem:

Theorem 4. For the double semion model in a honeycomb lattice Σ, the GE of spins for the ground state |+, + DS is given by
where n p is the number of plaquettes (faces) in Σ.
Proof. Remember the ground state |+, + for the toric code in the honeycomb lattice. This state can always be written as |+, + = c ijk... |ijk... with c ijk... ≥ 0 in the local basis {|+ , |− } for every spin. From here, it follows that | ijk.
where the equality holds if p i , q j , ... are all non-negative up to a common phase, and the maximum is achieved by choosing p i , q j , .. to be non-negative. Next, let us define d ijk... = c ijk.. e iφ , where φ is an arbitrary phase. Then it follows that | ijk.. p i q j ...d ijk.. | ≤ ijk.. |p i q j ....c ijk |. A maximisation over the left hand side gives the maximal overlap, say, for the double semion model by choosing appropriate phase factors and coefficients. But in any case, this will be smaller than maximising the right hand side whose maximum is, in this particular case, the maximal overlap for the toric code. From here the result just follows.
From the above Theorem it follows that the topological contribution is E γ = 1, which is again equal to the topological entropy for this model. Moreover, even if we do not discuss it here in detail, we expect the GE of blocks to obey similar laws as for the toric code model, since the ground state can also be disentangled by local moves [27]. Similarly, the entanglement properties of other ground states and excited states obtained by acting locally on |+, + , remain identical to those of |+, + (c.f. Lemma 1).

VI. COLOR CODE MODELS
Topological color codes [12] are defined as follows: qubits reside on vertices (not links!) of a lattice that is a 2-colex. A 2-colex Σ is a 2D lattice embedded in a torus of arbitrary genus g, with the following properties: (a) every vertex of the lattice has degree 3, i.e. precisely 3 edges meet at each vertex; (b) the faces of the lattice are 3-colorable, i.e. the faces can be coloured with 3colours (e.g. red, green and blue) in such a way that no two adjacent faces have the same color, see e.g. Fig 10. Two types of commuting operators are defined on each plaquette (or face) p, The Hamiltonian of the color code is thus For this model in 2D lattices, the ground space has degeneracy 2 4g , where g is the genus of the underlying Riemann surface. For instance, on the torus this degeneracy is 16. Let us define the group G X p as the group generated by the products of all possible plaquette operators B X p for all three colours. Notice that the product of all plaquette operators B X p corresponding to the same color c corresponds to the same element of the group G X p (namely, the action of a σ x operator on all the qubits). Therefore, we have the non-local constraint x . (32) The above relation implies that the size of the group is |G X p | = 2 np−2 , with n p the number of plaquettes in lattice Σ.
The ground level subspace of the color code model is constructed in a very similar way as we did for the toric code. Specifically, for the color code on a torus the ground space reads L = span{|i, j, k, l , i, j, k, l = 0, 1}, where Operators {w 1 , w 2 , w 3 , w 4 } above correspond to noncontractible loop operators, see Fig. 10.b. We do not enter here on the specific action of these operators on every spin in the loop, and just mention that it is local. The 16 vectors |i, j, k, l are orthonormal and stabilised by G X p , and therefore form a possible basis of the ground level space of the model on a torus.
As for the toric code, the action of loop operators corresponds to local unitary operations on a reference state |0, 0, 0, 0 , and therefore do not change the entanglement properties (c.f. Lemma 1). Therefore, from now on we shall consider only the entanglement properties of the ground state |0, 0, 0, 0 , defined as In what follows we show two different approaches to compute the GE for this model. The first approach is based, as all the calculations before, on finding upper and lower bounds for the GE. This approach works well for e.g. finding the GE of spin and blocks on regular lattices such as the honeycomb lattice. The second approach uses a different technique, building on recent work [28]. The main idea is to exploit that color code ground states are instances of so-called Calderbank-Shor-Steane states of self-orthogonal type. This refers to the fact that where S ⊆ Z n 2 is a collection of bit strings with a particular structure i.e. it is a self-orthogonal classical linear code. This property is in fact displayed by color code states associated with arbitrary lattices (2-colexes) and for surfaces of arbitrary genus. As a result, the second approach will lead to a general calculation of the GE for arbitrary color codes.

A. Bound method
Let us consider the color code on a honeycomb lattice Σ. For this lattice, and with this method, we can provide the following Theorem: For the color code model in a honeycomb lattice Σ, the GE of spins for the ground state |0, 0, 0, 0 is given by where n p is the number of plaquettes (faces) in Σ.
Proof. First, since the group G X p has size |G X p | = 2 np−2 , we obtain a lower bound on the maximal overlap Λ max ≥ |G X p | −1/2 , and hence an upper bound on the GE given by E G ≤ n p − 2. Second, in order to derive a good lower bound on the GE, we consider the natural bipartition for the honeycomb lattice, see Fig. 11.a. One can convince oneself easily that any group element in G X p must have nontrivial support contained in both the A and B sets, and therefore the subgroups G X p (A) and G X p (B) are both trivial and contain only the identity element. For the color code, the reduced density matrix of A also satisfies the property ρ 2 [29]. Therefore, Λ 2 max ≤ 1/|G X p |, and hence E G ≥ n p − 2. We therefore arrive at E G = n p − 2, which proves the theorem.
Some comments are in order. First, we see that the term of topological origin is E γ = 2, which is once again identical to the topological entanglement entropy [29]. Second, this topological geometric entanglement is twice that of the toric code, which is a direct consequence that the color code model we consider here is equivalent to two copies of the toric code under some given RG procedure [13]. Third, a similar result holds also for the ground state associate with the group G Z p generated by the product of B Z p operators. Fourth, excitations and other ground states related by local operations will have the same entanglement. And fifth, the GE of blocks will consist again of an area-law term plus a topological contribution which will be again E γ = 2. This can be seen, e.g. by choosing the blocking from Fig. 11.b, and appropriately disentangling qubits within each block.

B. Color codes and self-orthogonal classical codes
Here we provide a generalisation of theorem 5 which will allow to compute the GE of color code states in a much broader setting i.e. for arbitrary lattices (2-colexes) and arbitrary genus. The technique used to arrive at this generalisation is entirely different from the proof of theorem 5. It involves connections between color codes and classical coding theory which were found in the recent work [28].
We consider a color code defined on an arbitrary 2colex embedded in a surface of arbitrary genus. In anal-ogy with (34) we consider the state We will show here that the result in theorem 5 extends to all states |0 : Theorem 6. Consider the color code ground state |0 for an arbitrary 2-colex Σ embedded in a surface of arbitrary genus. The GE of spins for this state is Before proving the result, we develop some preliminary material. A linear subspace C ⊆ Z n 2 is called a (classical) binary linear code of length n. The elements of C are called its codewords. For every binary linear code C of length n, define an n-qubit state We shall refer to any state of this kind as a Calderbank-Shor-Steane (CSS) state. It easily follows from (37) that |0 is a CSS state. More precisely, we have |0 ≡ |C X p where the linear code C X p is defined as follows: writing every element g ∈ G X p as with u i ∈ Z 2 , the code C X p is simply given by the collection of all bit strings (u 1 , . . . , u n ) arising from elements of G X p in this way. In particular, one has which will be important below. It turns out that color code states are CSS states of a particular kind. A linear code C is called self-orthogonal if every two codewords are orthogonal: s T t = 0 for every s, t ∈ C. A CSS state for which the underlying classical code is self-orthogonal will be denoted CSS state of selforthogonal type. In [28] the following was shown: Lemma 3. Consider a color code associated with an arbitrary 2-colex Σ. Then the ground state |0 = |C X p is a CSS state of self-orthogonal type.
Next we provide a calculation of the GE of arbitrary CSS states of self-orthogonal type-by lemma 3 this will yield the GE of color code ground states. Note however that not all CSS states of self-orthogonal type are color code ground states i.e. our results extend beyond the color code setting. It is also interesting to point out that toric code ground states, even though they are CSS states, are not of self-orthogonal type.
Consider an n-qubit system and let |Φ = |φ [i] be an n-qubit product state, where Define an associated n-qubit tensor product operator as follows: Now consider an arbitrary self-orthogonal linear code C of length n and the associated n-qubit CSS state |C .
The following lemma, proved in [28], relates the overlap C|Φ to the expectation value C|A|C : Lemma 4. Let C be a binary self-orthogonal linear code of length n and let |Φ be an n-qubit product state. Then Lemma 4 will now be used to prove: Let |C be an arbitrary n-qubit CSS state of self-orthogonal type. Then the GE of spins is Proof. As before, let Λ max denote the maximal value of | C|Φ | when |Φ ranges over all product states. First, taking an arbitrary u ∈ C yields | C|u | = |C| − 1 2 so that we obtain the lower bound Λ max ≥ |C| − 1 2 . Second, since |C has real and nonnegative amplitudes in the standard basis, the maximal value of | C|Φ | will be reached by a product state which has real, nonnegative amplitudes as well. Consider an arbitrary product state |Φ with real and nonnegative amplitudes and the associated tensor product operator A as in (42). Using lemma 4 we find where · denotes the operator norm. Since each |φ [i] is a unit vector with real and nonnegative amplitudes, each A i is a real orthogonal matrix. Hence A i = 1. This proves the upper bound Λ max ≤ |C| − 1 2 which matches the previously found lower bound. The result now follows immediately.

VII. GENERAL FORMALISM FOR THE BOUND METHOD
In order to extend our results to the non-Abelian setting (e.g. quantum double models), it is worthwhile to generalise the bounding technique that we have used so far. Looking back, we have calculated the GE of quantum states with a certain (Abelian) group symmetry, and our bounding strategy involved two steps. First, we chose suitable product states to lower bound the maximal overlap (and hence to upper bound the GE). Second, we picked certain bipartitions and determined the maximal overlap with bipartite product states in order to upper bound the maximal overlap with multipartite product states (and hence to lower bound the GE). As we have seen, this procedure gives tight bounds in almost all the cases considered.
We shall now distill the essence of this bounding strategy and derive general bounds on the GE of quantum states with a non-Abelian group symmetry. Other than placing all of our previous examples under a common roof these bounds may even be of independent interest beyond the present study of topological order.

A. Overview
Let G be a finite group. Suppose its n-fold direct product G = G ×n acts on some m-partite where |Φ = |φ [1] ⊗ · · · ⊗ |φ [m] is some reference product state such that any two local states O t k |φ [k] and O u k |φ [k] are either identical or orthogonal.
We can define the global stabilizer G Φ ⊂ G as the subgroup which leaves the reference state |Φ invariant. Note that G Φ need not be normal in general, so the collection G := G/G Φ of cosets is not necessarily a group. It is clear that the state is correctly normalized. Picking any component O t |Φ we get a trivial lower bound on the maximal overlap Λ max (Ψ) and hence by (6) an upper bound E G (Ψ) ≤ n log 2 |G| − log 2 |G Φ |. Now let H = H A ⊗ H B a bipartition of the Hilbert space into subsystems A and B with the reference state |Φ = |φ A ⊗|φ B partitioned accordingly. For each operator O t we define a truncated operator O t X as the restriction of O t to H X and the identity in HX . We furthermore define the local stabiliser G Φ (X) ⊂ G for subsystem X as the subgroup whose truncated operators O t X leave the reference state |Φ invariant. It is clear that the global stabiliser is contained in any local stabiliser, so G Φ (X) := G Φ (X)/G Φ makes sense. Note that again G Φ (X) is not necessarily a group. Nevertheless, for bipartitions with the property we will find (see Lemma 6 below) that the largest eigenvalue of the reduced density operator ρ A equals |G Φ (A)| |G Φ (B)|/ |G| |G Φ | which by (7) immediately implies a lower bound on the GE. We can summarize both bounds on the GE in the following As is evident the derivation of this general bound is independent of the reference state |Φ or any physical model in which the state |Ψ might arise, including any geometry or topology that might be involved. So the problem of obtaining actual bounds for specific states in specific models with specific geometries (topologies) reduces to the much simpler problem of analysing the interplay between the group G, the global stabiliser G Φ and the local stabilisers G Φ (X) in the particular case at hand. Additionally, specific models will often allow us to characterise these stabilisers in purely geometric (combinatorial) terms, which provides a constructive way to find bipartitions obeying (48). Remarkably, this strategy even yields exact values for the geometric entanglement in many interesting cases as we will show (and have seen already). Namely, with a mild assumption in addition to (48) we immediately get the geometric entanglement of |Ψ is given by We would like to mention that if G is Abelian then both G and G Φ (X) are indeed groups and one can determine the spectrum of ρ A directly in terms of these. This has been discussed in Ref. [6] for the special case G = Z 2 . There our G is denoted by G and our G Φ (A) by G A (and similarly for the other subsystem B).

B. Some Details
Here we will supply some additional details glossed over before. First the precise definitions of global and local stabilisers.
This shows how we obtain arbitrary reduced density operators from the state (47).

Lemma 5 (Reduced density operator). For any bipartition H = H A ⊗ H B into subsystems A and B one has
Proof. Avoiding the sum over cosets, we may rewrite (47) as hence the global density operator reads By assumption {O t B |φ B | t ∈ G/G Φ (B)} is a set of mutually orthogonal states in subsystem B so we can use it to take the partial trace: Since the local expectation value is non-zero precisely for v ∈ G Φ (B) we obtain We turn to those bipartitions obeying (48).

Lemma 6 (Spectrum of reduced density operators). Let
Then the spectrum of ρ A is flat with non-zero eigenvalues Proof. From Lemma 5 we conclude that and now the observation is enough to prove the claim.

VIII. QUANTUM DOUBLE MODELS
Consider a planar graph Γ = (V, E, F ) with vertices V , edges E and faces F , and sizes |V | = n s , |E| = n and |F | = n p . The Hilbert space of the quantum double model is defined as H := H 1 ⊗ · · · ⊗ H |E| where each local H i CG. Vertex projectors (acting on stars) read where each individual vertex operator A g s acts on a vertex (star) s ∈ V by multiplying all edges meeting at s by g ∈ G from the left, provided all edges point towards s. Otherwise an edge is multiplied by g −1 from the right. It is clear that vertex operators on different vertices commute. Similarly, plaquette projectors are defined as which select configurations where the ordered product of group elements along an oriented circuit C p around p is the unit element of G. The Hamiltonian is the sum of these mutually commuting projectors over vertices and plaquettes, namely We identify the direct product group from our general discussion in the previous section as G := G ×|V | and obtain its natural action on the Hilbert space by collecting individual vertex operators into the joint vertex operator where t = (g 1 , . . . , g |V | ) ∈ G. Clearly, any joint vertex operator has product form: given a directed edge (s, s ) with value |x an element t = (. . . , g s , . . . , g s , . . . ) ∈ G is easily seen to act locally as We can obtain a ground state of the quantum double model by choosing the reference product state |e := E |e and projecting it on the common +1 eigenspace of all vertex operators In order to proceed we need to figure out the actual form of the global and local stabilisers in the quantum double model for our particular reference state. As for elements of the global stabiliser G e , it follows immediately from (56) that g s = g s iff vertices s and s are connected by an edge. In particular, we have G e G if the graph Γ is connected. Note that G e is not normal in G generally. As far as the local stabilisers G e (X) are concerned, we can equivalently define them combinatorially. Partition the vertices V into clusters V i (X) based on whether vertices are connected by an edge in E X . (That is, two vertices are in the same cluster iff they can be connected by a path which lies completely in E X .) Then t = (g 1 , . . . , g |V | ) ∈ G e (X) iff its components g s are constant on clusters. If E X connects all vertices into a single cluster then G e (X) G. It turns out that local and global stabilisers are related quite favourably: Proof. We want to cover the vertices V with a set {V i (X)} of overlapping clusters obtained from any of the subsystems A or B. Then if t = (g 1 , . . . , g |V | ) ∈ G e (A) ∩ G e (B) its components g s must be constant on each cluster, and hence constant on the whole of V because of the overlaps. Clearly, we can obtain a global cover by constructing a local cover for the vertices along each (edge) path γ.
So let γ be a path starting at vertex s and let V 1 (A) be the enveloping cluster of s as given by subsystem A. If γ never leaves this cluster we are done. Otherwise we may assume that s lies at the boundary of V 1 (A), hence we will reach a distinct cluster V 2 (A) by advancing a single edge (s, s ) along γ. Clearly, this edge must be in E B and thus there exists a cluster V 1 (B) containing both s and s . Then by induction {V 1 (A), V 1 (B), V 2 (A), . . . , V l (A)} is a local overlapping cover along γ.
So in order to find a suitable bipartition with G e (A) = G e (B) = G e we simply need to select a subset E A of edges such that both E A and E B = E \ E A have the single cluster property. Indeed, for the square, triangular and Kagome lattice such bipartitions exist, see Fig. 4(a) and Fig. 6 for examples. Using Theorem 9 we then obtain Theorem 10. For any quantum double model on a square, triangular or Kagome lattice Σ, the ground state |Ω has the geometric entanglement where α = log 2 |G|.
Note that |V | = n s (number of stars) is proportional to the volume of the surface.
From [24] and the proof of Theorem 2 it becomes clear that we can calculate the GE of blocks of linear size k ≥ 2 directly from the renormalized graphs Γ k shown in Fig. 2(e). On a square lattice of linear size 2k λ there are N k = 2λ 2 such blocks, and the renormalized graph Γ k has |V k | = N k (k + 1) vertices. We merely state the result for the GE of blocks.

Theorem 11.
For any quantum double model on a square lattice Σ, the ground state |Ω has the geometric entanglement of blocks of size k where α = log 2 |G|.
Clearly, |V k | is proportional to both the renormalized volume of the surface and the boundary area of a single block (with a natural value of 4k). This is exactly the generalisation of (14) to the non-Abelian case. As a final remark, we notice that the topological contribution to the geometric entanglement for quantum double models is given by E γ = log 2 |G|. This again coincides with the corresponding value of the topological entanglement entropy.

IX. Eγ AWAY FROM FIXED RG POINTS
So far we just discussed the topological contribution to the GE for specific ground states of models that turn out to be RG fixed points, and therefore representative of their respective topological phases. But, what if we are away from the RG fixed point? Is E γ robust under a perturbation? In this section we address briefly this question, and arrive to the conclusion that, indeed, E γ is a robust property of the topological phase.
There are several ways of checking the robustness of E γ under perturbations. One possibility is to perform a numerical analysis. In this sense, it should be possible to do large-scale calculations using tensor network methods. While we shall not carry out these calculations explicitly in this paper, we will explain different strategies based on tensor network algorithms in the Appendix. Nevertheless, we provide a small-size exact calculation at the end of this section. Another option is a perturbation theory analysis. Yet, a more intuitive alternative is the following argumentation based on RG fixed points: all the models considered here are RG fixed points and hence representatives of their respective topological phases. The longdistance properties of any state in one of these phases do not change under local RG transformations and, thus, are equivalent to those of the fixed point. This, in particular, is true for the long-range pattern of entanglement, and hence for E γ . Thus, any non-relevant and shortrange perturbation driving the Hamiltonian away from the fixed point will produce ground states with the same E γ . Nevertheless, we expect a change in the short-range pattern of entanglement, and hence in the bulk term corresponding to a boundary law. For short-range perturbations the change involves modifications of the pre-factor of the boundary law as well as the possible appearance of sub-leading O(L −ν ) corrections. Hence, Eq.(4) applies when away from the fixed point.
To double-check the above claim, we now perform a simple perturbation theory analysis of the robustness of the GE of the toric code model under external shortrange perturbations such as magnetic fields. This analysis provides upper and lower bounds for the GE, and complements the above argumentation on the robustness of E γ based on RG.
Let us then add a perturbation to the toric code Hamiltonian on the square lattice, and see how the ground state changes. For simplicity, we consider the case of an infinite plane, where the ground state |0, 0 is non-degenerate, and hence we can use non-degenerate perturbation theory. The perturbed Hamiltonian will be where H TC is the toric code Hamiltonian, V is the perturbation, and λ 1. Non-degenerate perturbation theory says that the new ground state can be approximated as where E 0,0 is the ground state energy and E φ,c is the energy of the excited state |φ, c, 0, 0 . We now consider the case in which the perturbation is an homogeneous magnetic field, e.g. in the x direction, (the case of z and y directions can be considered similarly). It is easy to check that in this case, the normalized perturbed ground state becomes with ∆ the energy gap to create a pair of flux and antiflux quasiparticles, and C = (1 + nλ 2 /∆ 2 ) −1/2 a normalisation constant.
Our aim now is to estimate the maximum overlap of the previous state with a product state of blocks of boundary size L. This can be done as follows: first, and as in the unperturbed case, we apply CNOT operations locally inside of the blocks so that qubits are disentangled in the unperturbed ground state. By doing this, we can focus on the entanglement of the state where |e k is the quantum state for the k-th disentangled qubit. The above equation is indeed equivalent to where S [j] x is the total spin in the x direction for the L spins in the boundary of block j, and x |e 1 ⊗ · · · ⊗ |e p .
Now we find upper and lower bounds to the maximum overlap of this state with a product state of the blocks. A lower bound can be easily obtained by the product state |0 ⊗(n−p) ⊗ |e 1 ⊗ · · · ⊗ |e p . Noticing that |e k is either |0 or |+ [24], and that the |+ contributions come only from qubits close to the boundary of the block, we have that for some positive ω = O(1) constant. The following up-per bound can also be found easily: Using the above bounds, one can check that for the GE we obtain, in the limit λ ∆ and L 1, The above equation is compatible with a leading change in the GE in the pre-factor of the boundary law. Also, the fact that both bounds leave the topological component E γ = 1 untouched seems to indicate that this is actually robust under the perturbation. Moreover, implementing finite-L corrections to these bounds provides O(L −ν ) corrections, which is consistent with our previous claims.
In order to further check the above arguments, we have computed numerically the GE for the toric code with n spins on the square lattice, exactly, for sizes up to n = 16, and where we added a perturbation that amounts to introducing a string tension in the Hamiltonian. More specifically, we perturbed the system by considering the PEPS tensors of the Toric Code ground state |0, 0 [37] and modifying the non-zero components that correspond to having "spin up" in all sites. That is, the components A i αβγδ of the perturbed-PEPS tensors are given by where the upper tensor index is the physical index and the rest are the bond indices (in the PEPS there is also another tensor like this but rotated 90 degrees, so that we have an ABAB... periodicity with a 2-site unit cell). For g = 0 one recovers the unperturbed toric code, whereas for g = ∞ one has a polarised state in the z-direction. This is the same kind of perturbation that was considered previously in Ref. [38], which adds a string tension with a somehow similar effect to adding a magnetic field in the z-direction to the Hamiltonian.
For such a perturbed topological state we did an exact calculation of the closest product state, for 1-site blocks (single spins) and also for 4-site blocks, up to n = 16 (i.e. up to 16 1-site blocks, and up to 4 4-site blocks). From this we estimated E γ by considering the GE as a function of the number of blocks, fitting it to a straight line, and extrapolating the number of blocks to zero. Calculations for larger systems required a significant amount of computational resources, and therefore could not be implemented with this approach (for these one would need e.g. the tensor network techniques explained in the Appendix). A summary of our results can be found in Fig. 12. Here we see that in the case of single-site blocks, the topological GE is not stable under the perturbation: it drops very quickly to zero without any sign of phase transition instead of staying close to E γ = 1 for a while. However, the 4-site block calculation shows that E γ also drops to zero, but quite slower than in the case of singlesite blocks. This is a clear indication that in the of 4-site blocks the topological GE is more robust under perturbation. We expect, thus, that as the size of the blocks becomes larger, the topological contribution to the GE becomes more robust. At this point we are constrained by the sizes in the exact calculations, but we take these results as a first-principle indication that the topological contribution to the GE tends to be robust under perturbation for sufficiently large block sizes, which is in fact compatible with the RG and perturbation theory arguments above.

X. CONCLUSIONS
In this paper, we have studied the GE in various topologically-ordered states that correspond to fixedpoints of RG. These are the toric code, double semion, color code, and quantum double models, in a variety of lattices. Generically, we found that the GE is typically composed of a boundary law term times the number of blocks plus a topological term. This is remarkable, since it is a signature of topological order in multipartite entanglement properties of the state, rather than bipartite. The topological term for the color code is twice that for the toric code, which is consistent with the recent result that the former model is equivalent to two copies of the latter [13]. Away from RG fixed-points, we argued that the same type of behaviour holds up to a possible subleading term, and thus the topological contribution is a robust property of topological phases of matter. The numerical evaluation of E γ for non-analytic cases is possible. This could be done e.g. in the context of Tensor Network methods (such as PEPS algorithms), as is usually done for the topological entropy [30]. We leave the specific implementation of these methods for future work.
In all the cases considered here, the topological GE turns out to be equivalent to the topological entropy. It would be good to prove, in general, whether this is always true or not, in order to understand if the GE can provide more information about topological order than the one that is already in the topological entropy. However, even if this were not the case, the GE would still be a useful tool to extract the quantum dimensions of the anyon model associated to a given topologically ordered state. For scenarios where reduced density matrices are hard to evaluate numerically, one may thus believe that looking at the GE may be more efficient. Moreover, the GE may also be a useful tool in order to extract minimally entangled states within the ground subspace, and hence the complete topological characterisation of the system [42].
Finally, we believe that other models could possibly be analysed with the techniques that we used in this pa-per. Specifically, similar results should also apply to the A-phase of Kitaev's honeycomb model [31] for which the toric code is (in some limiting cases) an effective model, as well as spin-liquid states with an emergent Z 2 gauge symmetry [32]. The possibility of studying the ground states of string-net models of Levin and Wen [5] in general, as well as topological quantum field theories [3], is also left for future investigation.
2.-Tensor-renormalization.-Another approach is to consider again blocks of size L×L, but where we compute a renormalized PEPS tensor for the block. This could be achieved by using different tensor-renormalization strategies, such as e.g. SRG [39] and HOSRG [40] but adapted to a finite 2d system with open boundary conditions. As a result of this tensor-renormalization, the topological quantum state is represented by a "renormalized" PEPS where each tensor corresponds to a block. Once this is computed successfully (which may be non-trivial), the necessary optimisation to compute the GE can be done using a product state directly for the renormalized single sites.