A Classical Model Correspondence for G-symmetric Random Tensor Networks

We consider the scaling of entanglement entropy in random Projected Entangled Pairs States (PEPS) with an internal symmetry given by a finite group G. We systematically demonstrate a correspondence between this entanglement entropy and the difference of free energies of a classical Ising model with an addition non-local term. This non-local term counts the number of domain walls in a particular configuration of the classical spin model. We argue that for that overwhelming majority of such states, this gives rise to an area law scaling with well-defined topological entanglement entropy. The topological entanglement entropy is shown to be log|G| for a simply connected region A and which manifests as a difference in the number of domain walls of ground state energies for the two spin models.


Introduction
For gapped models of local Hamiltonians on a finite dimensional lattice, the entanglement between a region A and the rest of the system is expected to be given by [3,4]: Here S(A) := − tr(ρ A log ρ A ) is the von Neumann entropy of the reduced density matrix ρ A of the groundstate in region A, ∂A is the boundary size of A, and α, γ ≥ 0 constants, both of which are independent of the choice of the region. The term o(1) stands for subleading corrections which go to zero with growing surface area of the region A. The first term in Eq. (1) shows that entanglement of the region scales with its boundary, instead of its volume as is the case for the majority of quantum states. It shows that the entanglement inherits the locality of the interactions, though we note that this is not always the case if the system is gapless. The second term γ is an additive correction to the area law scaling and is a signature of topological order in the system [11,12]. It is known as the topological entanglement entropy. For this quantity to be well defined, it is imperative that α is independent of the choice of the region.
However, Eq. (1) does not always hold. For example, we cannot have α independent of A for systems with disorder. For systems with translation invariance, it is possible that Eq. (1) applies generally. Even the problem of proving that S(A) ≤ O(|∂A|) is a major open question in dimensions larger than one (see [7,1,2] for the 1D case). Therefore, it is an interesting question to identify families of models for which one can argue Eq. (1) gives the right scaling of entanglement. One particular example are renormalization fixed-point models, such as string-net models [5,10,13]. There Eq. (1) holds with no correction term (denoted by o(1) there). However, it remains a challenge to go beyond such models. We note that numerically one can nonetheless verify the right form in many cases for the 2-Renyi entropy.
Given the difficulty of finding a general argument for the scaling of Eq. (1), one can ask whether it holds in the generic case, i.e. for most states in the class considered. If we consider ground states of gapped Hamiltonians as the class of states, this too seems like a very hard problem. Another approach is to consider states which automatically satisfy the weaker form of area law S(A) ≤ O(|∂A|) and ask whether generically the entanglement scales as Eq. (1). This is the path we take in this paper. In particular, we consider Projected Entangled Pairs States (PEPS), a useful class of tensor networks states with entanglement structure mimicking the locality of the interactions [17]. They are known to approximate well thermal states of local models in any dimension [8,14], and are expected to also approximate well the groundstates for gapped models (yet a proof for the latter is challenging). Seeing how the topological entanglement entropy manifests in such models may provide a clue towards proving the same in more general settings.
The case of random PEPS where each tensor is completely random was considered in an insightful paper by Hayden et al [9]. Their motivation was to study features of the AdS/CFT correspondence in a tensor network picture. Among other findings, they showed that the entanglement entropy scales as S(A) = log(D)|∂A RT | + o(1).
where D is the so-called bond dimension of the PEPS (the dimension of internal degrees of freedom of the state responsible for its entanglement) and |∂A RT | is the length of the Ryu-Takayanagi surface in the bulk for a boundary region A. The error term o(1) goes to zero as the bond dimension grows. The authors demonstrated a mapping from the problem of computing the entanglement entropy to the problem of computing the difference of free energies due to a domain wall in a certain classical Ising model associated to the state (where the log(D) becomes the inverse temperature of the model). Although the authors considered tensor networks in the relevant setting for holography (in which tensors on the bulk of the 2D space are fully contracted and only the tensors of the 1D boundary have physical degrees of freedom), their construction readily extends to the setting of PEPS, e.g. on a square lattice with one physical index on each tensor of the network. In order to extend this construction, we must acknowledge two major drawbacks of the original approach. The first is that the result relies on the "physical" boundary dimension being infinite. In the condensed matter setting, we care about systems of particles with finite degrees of freedom (e.g. qubits with physical dimension 2). This will require an additional constraint on our region A that its surface area to volume ratio be sufficiently small. Another drawback is that the topological entanglement entropy γ is always zero in the generic case of fully random PEPS. This is expected, as one expects to be able to adiabatically connect them to a trivial state. We work around this drawback by providing an additional constraint in the form of an internal symmetry on each tensor.
In more detail, in Ref. [16], Schuch, Cirac and Perez-Garcia gave a construction of PEPS with topological order associated to a finite group-incidentally a quantum double model-by imposing that the virtual degrees of freedom are symmetric under a regular unitary representation of the group. They also showed that for renormalization fixed-point PEPS, Eq. (1) holds with γ = log |G|, with |G| the order of the group associated to the topological order (e.g. G = Z 2 for the toric code).
In this paper we take an approach inspired by [9] and using the construction of [16] in an effort to calculate the form of entanglement entropy in random PEPS states with an internal symmetry associated to a finite group. We argue that generically, for certain regions A with surface area to volume ratio sufficiently small, regions which can always be found in the thermodynamic limit.
In analogy to Ref. [9], we also obtain our result by mapping the problem to computing the free energy of a classical model. Interestingly, in the case of non-trivial topological order, the model has (in addition to the usual local terms) a non-local energy term which counts the number of domain walls of a configuration. This is the term which gives the log |G| correction to entanglement at sufficiently low temperatures/high bond dimension.
The organization of the rest of the paper is as follows. In section 2 we review the construction of Hayden et al [9], noting important differences between our derivation and theirs. In section 3 we derive the classical model associated to a case with a local virtual symmetry, arriving at an Ising action with an additional non-local term. Then in section 4 we argue for the scaling of Eq. (3) analytically for the case of polynomial bond dimension (in the system size). Finally, in section 5 we summarize the results and provide some additional comments.

Second Rényi Entropy in a Random Tensor Network
In this section we review the construction of Hayden et al [9]. We choose to use mostly the same notation for clarity to readers familiar with this previous work. We will work through the calculation in necessary detail. This is for two reasons. First, we want to perform the calculation in a similar way though with several differences which are more convenient when the local virtual symmetry is introduced in the following section. Second, by expanding on the calculation in more detail, it becomes easier to comment on some key assumptions that are made.
We consider a two-dimensional square lattice throughout the entirety of this work for simplicity, though these results straightforwardly generalize to any lattice where the correspondence is then to a classical model on the associated interaction graph. The sites of the lattice are indexed by x. We associate a finite dimensional vector space to each of them. The tensor product runs over the neighbor sites y of each vertex x. To each neighbor we associate a D xy -dimensional vector space. Each site also has a d x -dimensional vector space V x,∂ , which will be associated with the physical degrees of freedom residing at each site of the lattice. Let D x be the total dimension of the space V x . Following the notation of [9], we also denote the indices related to physical degrees of freedom by so-called 'dangling' indices as they are not connected to other vertices.
Let |V x be a random state on V x , chosen from the Haar measure. Let |xy be a maximally entangled state on V y,v ⊗ V x,v , for two neighboring sites x, y. We start by considering a random (non-normalized) PEPS state defined as (see Figure 1): with the first tensor product running over all maximally entangled states acting on the vectors spaces V y,v ⊗ V x,v for all neighboring x, y, and the second over all sites. This leads to a generic form of the (non-normalized) density operator given by The subscript P indicates a trace over the vector spaces V x,v , corresponding to the contraction of internal indices. Note that ρ is a linear function of independent pure states living on the vertices. Our primary objective in this section is to compute the second Rényi entropy for a region A: where ρ A is the reduced density matrix of a region A.

Ising Model Correspondence for Topologically Trivial PEPS
Rather than starting from the entropy, we begin by stating the correspondence of [9] between the averaged trace over the choice of random tensor of a numerator and denominator (separately) of Eq. (7) and the partition function of a classical Ising model. This is presented in the following Lemma that will be a central component of this work.
Lemma 1 (Two-copy Trace/Ising Model Correspondence. From [9]). Let Z 1 := Tr ρ 2 A for a reduced density matrix ρ A on a region A which is described by a random PEPS. Then the Haar average of Z 1 equals the partition function of a classical Ising model of spins {s x } residing on the same lattice with Hamiltonian given by: where d x is the physical bond dimension of the tensor at x, D xy the virtual space bond dimension between tensors at x and y, and D x the total dimension of the tensor at x. The variable h x is the boundary pinning (magnetic) field which is 1 for x ∈ A and -1 for x ∈ A.
If we instead have an average over Z 0 = (Tr ρ A ) 2 , then we have the same correspondence but without a boundary pinning field, i.e. h x = 1, ∀x.
Proof. We begin by considering Z 1 (and Z 0 ) in a slightly more convenient form, where the second equality for Z 1 is simply a result of using the so-called 'swap trick'; the operator F A swaps the two copies of the system in region A only. Let's start with Z 1 . Insert eq. (6) into Z 1 and take an average to obtain, Note that the average of Z 1 follows through to just an average over the random vertex states |V x since the trace is a linear operator and the entangled states are not random. To compute this average, we take an arbitrary reference state |0 x so that our random states are just random unitary rotations from the reference, |V x = U x |0 x . We can then take the average by integration of unitaries with respect to the Haar measure [6], In this expression, I x denotes an identity operator and F x the swap operator each acting on the two copies of the Hilbert space at vertex x. F x is defined in a similar way as F A , swapping the two copies of the state at vertex x. Replacing the average in the expression with this result, we now have an expression consisting of 2 N terms for N vertices. Each of the terms indicates all combinations of either an identity or swap on each vertex. In order to simplify this expression, we introduce a classical spin variable, s x , where s x = +1 to indicate the use of an identity and s x = −1 to indicate the use of a swap operator. This representation forms the basis for our correspondence and leads us to consider Z 1 as a partition function of spins {s x } such that and The trace in eq. (14) is over an expression which is a tensor product of terms acting on different parts of the Hilbert space describing two copies of an N -vertex PEPS. The first part of the expression inside the trace consists of all maximally entangled density states describing internal lines in the network and acts trivially on the dangling indices. As for the swaps in the expression, F A acts only on the dangling indices in the region A and the remainder of the expression is either a swap or identity depending on the Ising variables describing each term. These swaps and identities act on physical and virtual spaces independently and thus can be considered as a tensor product of swaps/identities on each space.
Thus, we are motivated to split the trace as the expression inside is a tensor product of operators acting indepently on the Hilbert space of dangling indices and that of the virtual indices, The separate Hilbert spaces are indicated by the subscripts ∂ and v for dangling and virtual spaces, respectively. These subscripts are used in the same way for a swap operator on a vertex x, i.e.
Note that F x,∂ is a single swap operator acting on the dangling leg subspace of vertex x and F ⊗4 x,v is a tensor product of swap operators acting upon each virtual leg subspace associated to a vertex x.
The trace over dangling indices can be simplified by considering one operator which is either a swap or identity depending upon the combination of swaps in the expression. In order to do this, we introduce a variable h x , the same boundary field as in [9], The trace over the dangling Hilbert space is then over an operator which is a product of swaps depending whether x ∈ A or not and which term in the partition function we are concerned with. This allows us to write the trace as We use this particular expression to describe an operator which is a swap or identity as desired for a particular combination of s x and h x . The result is a trace over a product of operators acting independently on the dangling Hilbert spaces associated with each vertex, x. As such, the trace amounts to a product of traces on each space individually. Each of these traces is d 2 x if it is over an identity and d x if it is over a swap, where d x is the dimension of the dangling Hilbert space on the vertex x. This depends on s x and h x in the same way the operator does so that our trace above is simply This is the first ingredient in arriving at an evaluated expression for eq. (14) Next, we wish to compute the trace over internal legs. This trace is a product of traces over both copies of each edge in the tensor network, i.e. traces over the four virtual Hilbert spaces in which two copies of an internal leg reside. The swap operators here act independently on both copies of one virtual component of each tensor. Thus, we can rewrite this product of operators as a product over swaps applied on each edge in the tensor network, What we have done here is taken our original product of swaps and rewritten it as a product over edges instead of over vertices. We wish to maintain the dependence on the set of Ising variables, moving this instead to the exponent on each swap so that these operators are swaps in the instance s x = −1 (or s y = −1) and is an identity when s x = 1 (or s y = 1). Another way to think about this equivalence is by considering operators acting upon the clusters of Hilbert spaces surrounding each vertex and re-dividing them into sets of operators instead acting upon Hilbert spaces describing each edge in the tensor network.
Since the operators in this product of edges act independently on edge Hilbert space, we can reformulate the trace as a product of traces of operators on each internal space separately, The new subscript on the trace indicates that it is now a trace over an individual edge in the tensor network and not all virtual indices. The action of the swap operators is such that we are tracing over maximally entangled states connecting x to y and x to y (no swaps or two swaps) or x to y and x to y (one swap coming from s x = −1 or s y = −1). Here, the tick symbol indicates the second copy of the Hilbert space for a given subspace attached to that vertex. Thus, we can rewrite the tensor product of two swaps as a tensor product of one swap conditioned by s x and s y with an identity. If s x and s y are the same, then the operator is an identity and if they are different, then it is a swap: ) .
Here, we switch subscripts on the conditional swap operators to indicate more clearly which Hilbert spaces the operator acts on. We now have a trace over a maximally entangled state with a nontrivial operator being applied on one of the two Hilbert spaces on which the maximally entangled state is defined. We can then trace out the maximally entangled state, reducing the trace to a trace over just one virtual Hilbert space and its copy. This then resembles the dangling leg trace where we have either a swap or identity but now over a Hilbert space with different dimension, D xy , the dimension of the internal legs. Thus, Now, we insert eq. (22) into eq. (20) and combine this with eq. (18) to obtain the full expression with the traces computed, This then allows us to write to the Ising action as Z 0 follows exactly except where the swap F A is now an identity, thus making h x = 1 everywhere.
The correspondence constructed above provides us with the primary tool by which we can determine the form of the entropy. We elaborate in the proposition below. This and the lemma above are precisely the result of Hayden et al [9] that we will expand upon in the next section. Here, however, we have to make note of one essential difference. In [9], D xy = d x = D is taken for simplicity. In our case, we wish to consider systems with finite physical bond dimension corresponding to real particles. Thus, in general, our physical bond dimension will remain relatively small while our virtual bond dimension may be large. This places an additional constraint in the following discussion, requiring our surface area to volume ratio to be sufficiently small (though we can disregard this constraint if we also allow our physical dimension to grow).

2nd Rényi Entropy with Trivial Topological Term
We begin computation of the entropy by considering a connected boundary region A with reduced density matrix ρ A . The results below still hold for A consisting of disjoint regions though this only becomes significant when we introduce the topological constraint later on. The second Rényi entropy S 2 (A) is given by Eq. (7). We can consider the numerator and denominator separately and denote them by Z 1 and Z 0 , aligning with Lemma 1. The entropy average can be expanded in powers of fluctuations which has suppressed fluctuations at large enough bond dimension [9]. With high probability, we can then approximate the entropy with the averages performed separately over Z 0 and Z 1 The ' ' is an asymptotic equality as the bond dimension grows to infinity where the fluctuation term vanishes. We argue in a later section that this equality holds for polynomially scaling bond dimension. From Lemma 1, we see that Z 1 and Z 0 are both partition functions of classical Ising models so that the entropy is then a difference of free energies of these models, at an inverse temperature which is a function of bond dimension. At large bond dimension (low temperature), this difference of free energies can be approximated by a difference in ground state energies. In the Z 0 case, the ground state is simply a ferromagnet with all spins aligned in which the Ising action (Eq. 24) is a constant.
In the Z 1 case, the ground state is a bit more subtle. Normally, if we allow the physical bond dimension to take on a large enough value as with the virtual bond dimension (d x ≈ D xy ), the ground state is characterized by a flip of spins in the region A from the ferromagnetic field, introducing a log D xy |∂A| energy cost due to |∂A| non-aligned edges. However, in the case where we wish to keep d x small, we may have a region A where V log d x |∂A| log D xy in which it is preferential in energy to remain in the ferromagnetic state when flipping the magnetic field in the region A. This is because the energy cost due to non-aligned edges is far greater than a non-aligned spin in the region A. In order to account for this case, we then place the restriction that A be a region with sufficiently small surface area to volume ratio (number of sites in region A much larger than edges on the surface) so that the energy cost of non-aligned spins in the region A remains larger significantly larger than the energy cost of non-aligned edges around the surface area of the region. We expect that in such a case the Z 1 free energy will be as desired-flipped spins in the region A-below some critical temperature but we leave numerical evidence of this statement to future work.
So, with the small d x taken into account we can say the following: with high probability, the 2nd Rényi entropy of the region A scales asymptotically in bond dimension (inverse temperature) with the boundary of this region, where α is a constant dependent on bond dimension and A is a region of sufficiently small surface area to volume ratio. Note that if our bond dimensions D xy and d x are equivalent then we have an area law entropy in the large D limit regardless of the shape of the region A. Our computation here is performed somewhat differently from [9]; however, we likewise can obtain an expression for the entropy which is asympotically equal to a difference of free energies of the same Ising model without a bulk term. This bulk term can be introduced easily but only becomes relevant when considering a holographic tensor network. We will comment on this briefly towards the end of this work. Now, we move on to introduce an intrinsic topological symmetry that exists on the virtual indices in the tensor network and see what effect this has on our expression for the entropy.

Entanglement Entropy with G-Symmetry
In this section, our aim is to arrive at an expression for the Von Neumann entropy in the G-injective case by bounding it below by the second Rényi entropy and above by the Schmidt rank. We first work through the Rényi entropy following the previous section but with local symmetry and then present the final result below.
We consider a random tensor network, but constraining each vertex in the lattice so that the virtual indices remain invariant under operation by a unitary representation of a finite group G. We denote such representation by U (g) and require that it be a regular representation for reasons that will become apparent later in this section. It is worth noting that these representations must be such that the dimension of the representation or the total dimension of a a direct sum of such representations equals the bond dimension. The latter case may by of further interest though we do not consider it here. In any case, we will consider a tensor network-again residing on a square lattice for simplicity-made from states |ψ x on V x such that where the unitary representations, U (g), act upon the 4 virtual subspaces surrounding a particular vertex, the identity is acting upon the physical subspace, and g ∈ G are group elements. We represent this graphically in Figure 2. The more general case follows where we chose the symmetry such that a unitary representation lies on every virtual edge associated to a vertex.

Ising Model Correspondence for Topologically Nontrivial PEPS
With this in mind, we generalize our previous lemma to find a correction to the entropy corresponding to nontrivial topological order. The following proof will parallel that of the previous Lemma 1.
Lemma 2 (Projected Two-copy Trace/Modified Ising Model Correspondence). Let Z 1 = Tr ρ 2 A for a reduced density matrix on a region A which is described by a PEPS with local invariance given by a finite group, G (see Eq. 29). Then the Haar average of Z 1 equals the partition function of a modified classical Ising model of spins {s x } residing on the same lattice with Ising action given by where d x is the physical bond dimension of the tensor at x, D xy the virtual space bond dimension between tensors at x and y, and D x the total dimension of the tensor at x. The variable h x is the boundary pinning (magnetic) field which is 1 for x ∈ A and -1 for x ∈ A. The non-local term is characterized by the number of domain walls in a particular spin configuration, η, and the cardinality of G.
If we instead have an average over Z 0 = (Tr ρ A ) 2 , then we have the same correspondence but without a boundary pinning field, i.e. h x = 1, ∀x.
Proof. To begin, we start with the same random tensor network construction as before but attempt to constrain our quantum state to one where local tensors are G-symmetric. We impose the desired symmetry constraint by restricting the random tensors to symmetric subspaces. To do this we first introduce our symmetric subspace projector, We can then project the states living on each vertex onto the symmetric subspace using this projector to obtain a slightly modified expression for the numerator of the entropy, Z 1 , from the previous calculation. We will do this on a uniform square lattice (n x = 4, ∀x) for simplicity, but the formalism extends to any random PEPS tensor network in terms of operators acting upon virtual subspaces in each edge.
Here we do not have the adjoint of the projector on the right side of the average as it commutes through the average. This average is done in the same way as before, leaving us again with an expression which no longer contains randomized vectors, From here our we wish to make the same attempt at a calculation we did in the previous section with the hopes that the added projector will not introduce many difficulties. First, we'll introduce the same Ising variable as before which allows us to write where eq. (13) still applies. Since P x acts only upon the virtual indices of a tensor,we can write the expression inside the trace again as a tensor product of operators acting independently on the virtual and physical Hilbert spaces, This first trace is identical to the previous case, so eq. (18) still applies. The second trace, however, is complicated by additional terms coming from the symmetric subspace projector. In order to break apart the second trace as a trace over individual edges in the tensor network, we first need to look at the projector we have introduced. This operator consists of unitaries which are independently applied to the virtual Hilbert spaces surrounding a particular vertex. As such, we can split the tensor product of all such operators into operators acting upon both copies of a given virtual Hilbert space. Thus, Essentially, what we have done is written our product of projectors and conditional swaps as a product of operators acting upon Hilbert spaces associated with a particular edge in the tensor network and its copy. This includes the unitary representation on each virtual Hilbert spaces as well as the swaps applied. We have rewritten the sum as a sum over four independent group elements chosen from G to specify that we have all terms where each virtual subspace in an edge and its copy may have a different group element attached to it. The N in the prefactor is simply the number of vertices in our tensor network, coming from a subspace projector being applied to each vertex.
We can write this expression a little more compactly in terms of the trace we wish to compute, leaving where the argument of the trace we wish to compute is The conditional swap operators here work in the same way as before, either flipping the spaces connected by maximally entangled states or not. We can thus combine this into one swap operator acting on either the xx subspaces or the yy subspaces conditioned by s x and s y , This then leads us to two possibilities. If no swap is applied, then maximally entangled states reside in the x and y virtual Hilbert spaces and x and y Hilbert spaces separately. Thus, our expression inside the trace becomes simply a tensor product of group elements U (g x g −1 y ) and U (g x g −1 y ) acting on disjoint Hilbert spaces after tracing out one virtual Hilbert space and its copy. The inverse comes from moving the operator from one Hilbert space to another through the maximally entangled state and multiplying group elements together.
On the other hand, applying the swap connects maximally entangled states from x to y and y to x . Now, when we take a partial trace we are left with just one operator acting on a virtual Hilbert space and its copy but not independently. This operator is given by U (g x g −1 y g x g −1 y ). We can consolidate these two possiblities in the follow equation, The traces are evaluated as Kronecker deltas because we have a regular representation (or direct sum of regular representations if the bond dimension is large) such that the trace over any non-identity element is zero. This introduces a subtlety into the available states since the trace for s x = s y can be nonzero for some nontrivial group elements residing on the lattice. Note that we still obtain the same factors D xy or D 2 xy as before. Returning to the original expression, we now have a product of delta functions dependent on our Ising variables which is being summed over group elements corresponding to each vertex in the tensor network. In order to fully evaluate the expression, we then need to count the number of terms in the sum over group elements for which the product of delta functions depending on the Ising variables is nonzero.
In order to perform this counting, we first partition {s x } into subsets given by the number of domain walls, η, in the lattice. This allows us to determine the number of factors in the product.
Let's first consider if we have one cluster (all s x = 1 or all s x = −1) then no swaps are applied anywhere and we only have conditions of type δ xy δ x y . These conditions tell us that vertices connected by edges on the tensor network must have the same group element associated with it, g x = g y for a given link. Thus, the only term which are nonzero in the sum are the ones for which all vertices in the tensor network have the same group element applied to each virtual index. This is independently and similarly true from the second copy of the tensor network. Thus, we have |G| 2 total nonzero terms in the sum since we have |G| possible group elements residing on each vertex in the tensor network as well as its copy. Now, suppose we add one domain wall to the Ising lattice describing out tensor network so that we have two clusters. This changes terms associated with edge crossing the domain wall to four-index deltas, δ xy x y . These delta functions are associated with link connecting two regions in which the lattice is otherwise homogenous in group elements. Each of these δ's associated with the same domain wall describe the same condition, g x g −1 y = g x g −1 y . This condition is the same for each edge crossing the domain wall since group elements outside the domain wall are the same (as well as inside the domain wall). Unlike the other type of factors in the product, this one only restricts one of the four groups elements associated with a link. This then gives us the freedom to determine one more set of group elements within the tensor network, adding a multiplicative factor of |G| to the number of nonzero terms. Thus, if we have η domain walls, then there are |G| 2+η nonzero terms in the sum. So, returning to eq. (37), Now, combining this expression with eq. (18), we obtain This gives us the Ising action We can write this in a slightly more useful form in which we can collect several terms as a constant that will be the same in both free energies associated with Z 1 and Z 0 , As before, we get the same expression for Z 0 but with a h x = 1 everywhere. This is the Ising action for our locally symmetric random tensor network and forms the basis of a modified Ising model which lends itself to computing the Rényi entropy which we do below.

Entanglement Entropy with Nontrivial Topological Term
As in the proof of Proposition 2, in the limit of large bond dimension, the entropy is approximated by the difference in free energies of two different (now modified) Ising models, We again approximate this difference of free energies as a difference of ground state energies which we expect to be the case above a critical bond dimension. Now, however, we will be apply Lemma 3 to determine the ground state energies of the corresponding model. In the Z 0 case we still have a ferromagnetic phase where η = 0 so that the ground state energy remains a constant. Z 1 , however, has flipped spins in the region A again assuming the surface area to volume ratio is sufficiently small in the case where we don't allow d x to be large. We expect this to be the case but again leave numerical evidence to future work. If A is composed of m disjoint regions, then the ground state contains m clusters of flipped spins such that η = m. Thus, the ground state energy scales with the boundary of the region A as before but now contains an additional term given by −m log |G|. The entropy thus scales asymptotically as the boundary with a negative logarithmic correction where α is again a function of the bond dimensions, Note that in the case that A is one simply connected region then the correction is − log |G|. It is also worth noting that this theorem still relies on the shape of the region A in the instance where our bond dimensions differ. Nonetheless, we can see how the topological entanglement entropy manifests in the entropy average as a difference in domain walls in the classical model ground states.
We can now present an expression for the Von Neumann entropy by way of an upper and lower bound. The Von Neumann entropy is lower bounded by the second Rényi entropy so that it follows immediately from Proposition 4 that Next, we wish to upper bound the Von Neumann entropy. We do this by first bounding the maximal Schmidt rank. Consider a state |ψ described by an arbitrary G-injective PEPS tensor network. We then consider a bipartition of the physical Hilbert space into regions A and A. We can then act on the region A by a linear operator M A composed of local operations which purifies the degrees of freedom in A. This produces a state |φ = (M A ⊗ I) |ψ whose Schmidt rank is greater than or equal to that of ψ state since M A is acting only on region A by local operations and can therefore only increase the maximal Schmidt rank. However, since the state is G-injective, we have an invariance given by this bipartition, |G| −1 g∈G U (g) ⊗|∂A| |φ = P ∂A (G) |φ = |φ . Our maximal Schmidt rank is then upper bounded by the rank of this projector S(A) ≤ max (log (rank ρ A )) ≤ log (rank P ∂A (G) |φ ) = log (D |∂A| xy /|G|) = |∂A| log D xy − log |G|.
(48) Thus, we have shown that our Von Neumann entropy is upper bounded by a general area law form in the G-symmetric case and is similarly lower bounded, though asymptotically so. Our Von Neumann entropy then asymptotically has the desired tight bound, where α is a constant dependent on bond dimensions, and |∂A| is the size of the boundary of region A.

Bound on Fluctuations at Finite Bond Dimension
In this section we attempt to bound fluctuations by reviewing the analysis done by Hayden et al [9], making a detour to employ necessary physical assumptions that restrict the bond dimension scaling by using a low-temperature expansion of a partition function. In order to simplify the problem, we assume the Ising model in question sits on a square lattice. The result, however, will hold in general for a model on any finite graph. We also take D xy = d x = D for additional simplicity but expect a similar result to hold if we keep d x small. We make this precise in the following Theorem. Proof. We will go about bounding this difference with some probability by way of Markov's inequality. In this proof, we inherit the previously mentioned definitions of the partition functions and Ising actions associated so the entropy. In addition, we denote the Ising action of the minimal energy configuration with boundary field h x by A 1 ] to be the large D limit of the partition function, we can write the exact nth Rényi entropy as Before employing Markov's inequality, we first choose to bound fluctuations away from the finite D partition functions by first considering From here we make use of the same insight given by [9] which is that the second moment, (Z (n) 1 ) 2 , can be interpreted as the partition function of a Sym 2n -spin model. The boundary field, h x = (1...n)(n + 1...2n), consists of two separate n-cycles in the region A and is an identity elsewhere. For large bond dimension, the minimal energy configuration is the same but this higher order model has energy where the constants are given from eq. 67 and are a difference of each normalization, hence the term in the above equation. We then continue from eq. 52 so that This recovers the result of previous work and is where we will continue with a somewhat different approach.
We continue by expanding (Z (n) for low temperatures. Since changing spins increases the number of domain walls, with A a more tractable Ising action given by From here we can expand the exponential by considering low energy configurations of the partition function. This is where having a square lattice greatly simplifies the problem. The lowest excited states are those with just one spin changed from the ground state, of which there are 2n − 1 values it may take and N sites that could be changed. Note that changing a spin next to the boundary of region A incurs a slightly different energy cost from changing a different spin, but in any case the energy incurred is proportional to log D. We can similarly consider "dimer" terms where we change adjacent spins in the same way and then even higher energy configurations with two non-adjacent spins changed or adjacent spins changed differently. This allows us to write where a i and b i are small constants with a 0 = 1 and b 0 = 0. This allows us to continue from eq. 54 such that (Z We want this quantity to be arbitrarily small for the Rényi entropy to be approximated within arbitrary precision, which is the case if D >> nN , i.e. if the bond dimension scales with the volume of the system (as well as the order of Rényi entropy). If we write the sum as an exponential, we have (Z where c n is a constant proportional to n. Then, by Markov's inequality, we write where we get an equivalent expression for Z Finally, we can bound the deviation from the exact area law by where in the second to last inequality we make use of the fact that δ ≤ 2 so that log (1 ± δ/4) ≤ δ/2. This bound holds with probability at least 1 − e Dc/D with D c = log (32/δ 2 )c n V .
Thus, we have demonstrated that our formula for the entanglement entropy holds to arbitrary precision for a bond dimension scaling (polynomially) with the volume of the system. We can therefore always choose our bond dimension to be much larger than the volume of the system so that the leading order term grows logarithically with the volume of the system. This implies that the area law effect is generically apparent for regions A which are volume logarithmic in size. We note also that the constant correction arising from symmetry constraints plays no role in this bound since its value does not change with bond dimension.

Conclusions and Future Work
In this paper we extended the method of [9], providing a connection between the computation of Rényi entropy with a difference of free energies of an Ising model. We generalize this method to include a local symmetry at each tensor in the tensor network lending itself to a non-local term in the Ising model and correspondingly a topological entanglement entropy present in the area law. We argue that this form holds exactly at zero temperature as well as small, nonzero temperatures scaling polynomially with bond dimension. While a specific case is presented, we expect these results generalize to any dimension and for a symmetry given by any local, finite group. We can also see how the topology of region A presents itself in the topological term, adding a power to |G| in the logarithm by enforcing more domain walls in the nontrivial free energy expression. Thus, we have shown that a random PEPS tensor network subject to a local symmetry has an area law with high probability at sufficiently large bond dimension (low temperature). This presents a larger class of states known to have area law than before as we have extended beyond mere renormalization fixed-point PEPS.
In future work, we would like to understand more precisely the phase transition in the classical spin model we consider. This will tell us precisely the bond dimension above which we may choose a random PEPS with area law with high probability. This analysis can be done by performing extensive numerical work on the classical, modified Ising model presented above. We note that this is computationally challenging due to the nonlocality of the modification, though not prohibitively so. We would also like to extend our model to a larger class by including more symmetries. For example, one may encompass all string-net models by way of generalizing the symmetry to include MPOs symmetries [15]. This would require determination of a symmetric subspace projector for any MPO symmetry or another way to represent a generalized symmetry on an otherwise random tensor network. Another avenue of interest may be to extend the G-symmetric random tensor network construction to one where additional tensors without associated physical Hilbert spaces are added to the network. This may be of particular interest in the case of understanding G-symmetry in a holographic tensor network. It also appears straightforward to extend the correspondence to Projected Entangled Simplex States whereby our entanglement occurs on simplices constructed from multiple vertices rather than between adjacent vertices.

A Higher Rényi Entropies
In this section, we will extend our results to a computing higher Rényi entropies, again paralleling results from [9]. This is mostly done to demonstrate how the topological correction manifests in higher entropy expressions. We will again use similar variables for clarity and will indicate a few differences as they appear.
We generalize the second Rényi entropy to the n-th Rényi entropy defined by We choose variables Z

Classical Model Correspondence for Higher Rényi Entropies
We can write a generalization of Lemma 3, extending the correspondence of these traces to partition functions of classical spin models. The proof for this Lemma will follow the same structure as for previous Lemmas but will make use of some examples to elaborate on the analysis. = Tr ρ n A for a reduced density matrix on a region A which is described by a PEPS with local invariance given by a finite group, G (eq. 29). Then the Haar average of Z (n) 1 is equivalently described by the partition function of a modified classical n-spin model of spins residing on the same lattice with Ising action given by where d x is the physical bond dimension of the tensor at x and D xy the virtual space bond dimension between tensors at x and y. The function χ counts the number of cycles of a permutation group element. Γ x is the permutation group element corresponding to the value of the spin on the same vertex. The variable h x is the boundary pinning (magnetic) field which is an identity on the n copies, I x , for x ∈ A and a cyclic permutation on the n copies, C (n) x , for x ∈ A. The non-local term is characterized by a permutation group element associated to a particular domain wall, η i , between clusters, Γ η i , and the cardinality of G.
If we instead have an average over Z (n) 0 = (Tr ρ A ) n , then we have the same correspondence but without a boundary pinning field, i.e. h x = I x , ∀x.
Proof. We first consider Z (n) 1 which can be modified by a generalization of the "swap trick" where we instead permute the n copies of the state cyclically in the region A, the action of which is given by the operator C (n) A . It's important to remember that the action of this variable is cyclic which will be relevant later in this section. With this in mind, we can write Our average on Z 1 is now an average over n copies of |V x V x |. This average is a symmetric subspace projector of the n-fold tensor power Hilbert space [6] so that where we replace the g x notation from previous work with Γ x to specify group elements coming from the n-copy average and not elements of our local symmetries groups. This is to avoid confusion between the two and keep consistency with our former definition for g x . These Γ x are permutation group elements acting upon n copies of a single site Hilbert space but can be non-cyclic group elements unlike the single site operation by C (n) A . The normalization constant remains the same as in previous work, representing a sum of traces over each group element in Sym n acting upon each single site Hilbert space of vertex x. From here, we can split the trace into virtual and physical Hilbert space traces as done in Sections 3 and A since the permutations still act independently on these Hilbert spaces and the local symmetric projectors are still acting on the virtual Hilbert spaces. The physical trace then shows up in the Ising action in the same way as the previous result [9]. We will give the expression for this term precisely when we present the full expression for the Ising action below.
Before we do this we must first deal with the trace over virtual Hilbert spaces. The permutation group elements act in a similar way to the swap operators in the analysis in section 4. Rather than expressing them as operators depending on Ising variables, however, we leave the operators in the general from of Γ x to indicate that they are some particular permutation group element acting upon all copies of a particular virtual subspace of x. This allows us to rewrite equation 37 as where Γ ⊗4 x,v represents the action of Γ x of the four n-copy virtual Hilbert spaces at vertex x. Note the additional sum which runs over all combinations of Γ x and Γ y . Since we are now working with n-copies of our state, we have a local symmetry projector acting on each of the copies resulting in a prefactor of 1 |G| nN and we have n group elements associate with each copy for both sides of a particular edge on the tensor network. The Ω in the above expression is constructed in the same way and is given by Performing a permutation on the virtual Hilbert spaces copies of vertex x and on vertex y is the same as performing a different permutation on x or y and an identity on the other. This new permutation on x is given by a composition of Γ x,v and Γ −1 y,v , the inverse group element. Since this is still just an arbitrary group element we can just leave Γ x,v and make Γ y,v an identity, mirroring eq. 39.
Since composition of group elements is again a group element, we can reduce the sum over Γ x and Γ y to just a sum over all possibilites for Γ x where Ω becomes Ω = |xy xy| ⊗n (U (g x 1 ) ⊗ · · · ⊗ U (g xn ) ⊗ U (g y 1 ) ⊗ · · · ⊗ U (g yn ))(Γ x,v ⊗ I). (70) We can now consider different cases for different permutation group elements, Γ x,v . The action of this operator permutes the group elements in the virtual Hilbert space and its copies associated with one vertex of a given edge. The reduction given by the above equation is then given by considering the permutation of group elements on each side as a relative shift in the two given by a single permutation of group elements on one side. Since the maximally entangled states associated with each edge are attached to a particular copy of the edge Hilbert space, permutation of operators on these spaces then allows us to push the unitaries on one side of the edge through as inverses in every possible configuration.
At this point we introduce what is called an edge group element, e i = g y −1 i g x i , to make our expression simpler. It is easy to see why this is useful if the permutation is an identity since we're doing no more than pushing together each copy of the edge independently. If multiple such copies of an edge are dependent, then such that we get a unitary of a product of these edge group elements. For example, if the permutation is a flip on two copies of a Hilbert space and an identity on the rest then we can write some unitary operator U (g x i g −1 y j g x j g −1 y i ) for some i and j tensored with a unitary over one edge element for every other copy. Since a trace is invariant under cyclic permutation of operators within the trace, it follows that the nontrivial unitary can also be written as U (g −1 . This can be extended to any number of dependent copies of an edge such that we can always write our unitaries in this way. Note that these unitaries are unique up to cyclic permutation of edge group elements. With this in mind, Γ x,v then connects copies of Hilbert spaces on one side of an edge to copies on the other side which gives rise to tensors of all possible unitaries over edge group elements where each edge group elements shows up only once in the expression. This allows us to define a map from Tr <xy> [Ω] to Tr x [Ω ] given by translating the cycles associated with a particular permutation group element to unitaries over edge group elements corresponding to these cycles. Here, the trace over x indicates a trace over the Hilbert space and its copies for one side of a particular edge (we've brought operators through maximally entangled states and reduced the other side to an identity). As a specific example, suppose Γ x,v is given by the cycles (23)(145) for the group Sym 5 where (ijk) describes a cyclic permutation of operators in Hilbert spaces i, j, and k (in this case we have ordered an edge Hilbert space and its copies in some way). Then the instance of Ω containing this particular permutation can also be written as Ω = U (e 2 e 3 ) ⊗ U (e 1 e 4 e 5 ) where each unitary operator acts on the Hilbert space associated with its particular edge group elements. Note that these unitary operators are unique up to cyclic permutations just as the indices in the cycles of a permutation group element are.
Enumerating these possibilites for large order permutation groups is burdensome, but we present the cases for Sym 4 as an example. Here, i, j, k, l ∈ {1, 2, 3, 4} with none of these indices being equal.
We see that considering the trace this way again gives us a product of Kronecker deltas as before. Returning to the connection to an Ising model, we have the same Sym n -spin model with the same boundary pinning field as in (Hayden et. al.) but we now obtain additional factors of |G| depending on solutions to equations given by the product of delta functions. It is important to note that the expression is otherwise identical to the previously obtained result. The higher order entropy permutation affects the number of nonzero terms in a sum of expressions that are otherwise unchanged by the topological symmetry.
There is a subtle difference, however, in that the number of additional solutions across each domain wall in the Sym n -spin model is not simply |G| but a multiple thereof. Nontrivial unitaries in the trace (those containing more than one edge group element) give rise to Kronecker deltas with additional solutions. Each nontrivial Kronecker delta gives another set of |G| solutions so that the additional factors of |G| introduced at each domain wall depends on the relative permutation across the domain wall. We quantify the number of these extra solutions by using χ(Γ η i ) to indicate the number of cycles in the relative permutation across a domain wall indexed by η i .
The number of added solutions at a given domain wall is then n − χ(Γ η i ). This coincides with no additional solutions if there is not a domain wall (χ(I) = n), one additional solution if the relative permutation is a swap on two copies (χ(F ij ⊗ I) = n − 1), and n − 1 additional solutions if the relative permutation is a cyclic permutation (χ((1 · · · n)) = 1). The prefactor can then be written in a slightly more complicated form as before so that The product over η i is a product over all domain walls in a particular configuration, {Γ x }, of the Ising model lattice. Notice how the separate product over bond dimensions translates to the same result as previous work, but we now have a power of |G| translating to the topological entanglement entropy as we will we see more clearly below. In a more concise way, we are considering a particular configuration, determining the relative permutation across a domain wall, and then constructing an equivalence to a system of equations given by a product of Kronecker deltas. We can systematically count the number of solutions to determine the addition factors of |G|. From here we can write the Ising action as with the boundary pinning field, h x , given by x , x ∈ A. (74) Grouping the constant terms in the Ising action, we recover eq. 64. For Z (n) 0 , the same analysis holds with a trivial boundary field whereby h x is an identity everyhwere.

Higher Rényi Entropies for a Nontrivial Topology
Now, we attempt to write down the form of a generalized Rényi entropy of a random PEPS subject to a non-trivial topology. Applying Lemma 5, we can again write the entropy as a difference of free energies, where the superscript on the free energies means that we are referring to the free energy of our n-spin model. At large bond dimension, or low temperature, the free energies are approximated by ground state energies. To determine the ground state, we note that the internal bond dimension term in the Ising action is zero as Γ x = Γ y , encouraging a ferromagnetic state in which the 2n-spins are aligned. The boundary terms vanish when the particular cycle is aligned with the boundary magnetic field so that Γ x is a cyclic permutation within the magnetic field and I otherwise. This encourages a uniform identity state without a magnetic field and exactly m domain walls around the region A with a magnetic field. The relative permutation across these domain wall is a cyclic permutation so that the sum over domain walls is simply 1 − n. Thus, Without the magnetic field, there is a ferromagnetic phase in which the ground state energy is a constant. So, we can write the difference of free energies and therefore the entropy as This is precisely the area law form as desired. We do this in the case where d x = D xy = D but expect a similar result to hold if we wish to keep d x small where we again require a sufficiently small surface area to volume ratio.