Superselection-Resolved Entanglement in Lattice Gauge Theories: A Tensor Network Approach

Lattice gauge theories (LGT) play a central role in modern physics, providing insights into high-energy physics, condensed matter physics, and quantum computation. Due to the nontrivial structure of the Hilbert space of LGT systems, entanglement in such systems is tricky to define. However, when one limits themselves to superselection-resolved entanglement, that is, entanglement corresponding to specific gauge symmetry sectors (commonly denoted as superselection sectors), this problem disappears, and the entanglement becomes well-defined. The study of superselection-resolved entanglement is interesting in LGT for an additional reason: when the gauge symmetry is strictly obeyed, superselection-resolved entanglement becomes the only distillable contribution to the entanglement. In our work, we study the behavior of superselection-resolved entanglement in LGT systems. We employ a tensor network construction for gauge-invariant systems as defined by Zohar and Burrello (2016) and find that, in a vast range of cases, the leading term in superselection-resolved entanglement depends on the number of corners in the partition, that is, corner-law entanglement. To our knowledge, this is the first case of such a corner-law being observed in any lattice system.


I. INTRODUCTION
Gauge theories are central in fundamental physics, originated in continuous field theory.In a gauge theory, the discussed system is subject to a local (spacetime dependent) symmetry, which give rise to new fields referred to as the gauge fields.In lattice gauge theory (LGT), [2][3][4] the continuous gauge theory is limited to a discretized lattice space (or spacetime).
The entanglement in LGT models is tricky to define, due to the nontriviality of the Hilbert space [31][32][33][34][35][36][37][38][39][40][41].However, when one restricts themselves to a single symmetry sector, also denoted by a superselection sector in LGTs, the discussion of entanglement becomes natural, as we explain in Sec.IV.In this paper, we study symmetry resolved entanglement [42][43][44] which have been extensively studied in many systems and led to interesting discoveries .The local symmetry imposes a large number of symmetry sectors, referred to as superselection sectors.The study of gauge symmetry resolved (SR) entanglement is specifically interesting, as it allows to separate the entanglement into two contributions: the SR entanglement, which is the only distillable entanglement when the gauge symmetry is fundamental, hence has to be obeyed by any operation performed on each subsystem; and the entanglement stemming from the division into different sectors, which detects topological effects.
We use the special structure of gauge-invariant states on TN to obtain a corner-dependent term in the superselection-resolved entanglement and discuss the various cases in which this is the leading term, that is, the entanglement's behavior follows a corner-law.Since studied systems typically have a constant number of corners (e.g., rectangular system with four corners or a strip on a cylinder with zero corners), the superselection-resolved entanglement in these cases becomes independent of system size.When the harvest of entanglement is desired as a resource for quantum technology, one may use the corner-law to find the partition with maximal entanglement based on the number of corners it contains.The rest of this paper is organized as follows: In Sec.II, we present the basics of LGTs, while focusing on pure gauge models on a square lattice.In Sec.III, we present the TN construction we work with, projected entangled pair states (PEPS), discuss its relation to entanglement, and present the PEPS construction for gauge invariant states on which our work relies on as proposed by Zohar and Burrello [1].In Sec.IV, we discuss symmetryresolved entanglement and its meaning in the presence of a local gauge symmetry.The reader who is familiar with PEPS and symmetry-resolved entanglement may jump straight to Sec.IV A, in which we bring all of the above together and obtain our main result regarding the behavior of superselection-resolved entanglement in PEPS-representable gauge-invariant states with an Abelian gauge group; the non-Abelian case is discussed in Sec.IV B. In Sec.V, we present numerical results on pure Z 2 gauge models, which are in line with our results and perhaps raise questions regarding the relation of entanglement and confinement in gauge-invariant states.
An important notation clarification is required before proceeding with the paper: Note that in the standard notation of entanglement of two-dimensional systems, arealaw (volume-law) refers the entanglement polynomial de-pendence on the boundary (bulk) size of the system, while in the standard notation of confinement (as explained in Sec.II), area-law (perimeter-law) refers to a decay that depends on the bulk (boundary) size of a Wilson loop.In this paper, we will use both notations according to the discussed property and also add the bulk/boundary notation, for clarity.

II. GAUGE-INVARIANT STATES ON A LATTICE
Gauge theories were first introduced in the context of quantum field theory: A local (i.e., time-and spacedependent) symmetry gives rise to so-called force-carrier fields or gauge fields, which mediate the interaction between matter fields [77][78][79].
LGT was later introduced as a numerical discretization tool for the study of continuous models [2][3][4].However, LGTs turned out to exhibit interesting properties on their own, and are now useful in the study of topological effects in many-body physics and as candidates for quantum error correction codes.
As in the continuous gauge models, LGT models consist of matter degrees of freedom (DoF) and gauge DoF.The matter DoF reside on the lattice sites, and the gauge DoF, mediating the interaction between two lattice sites, reside on the edges between the lattice sites, as demonstrated in Fig. 1a.In this paper, we focus on twodimensional square lattices, but the generalization could be done to higher dimensions, nonabelian symmetries, or different lattice types.We follow the Hamiltonian formalism of LGT in Ref. [3], i.e., only space is discretized and time remains continuous.G denotes the group corresponding to the local symmetry.For simplicity, we focus here on Abelian symmetry groups (for the generalization to nonabelian symmetries see, e.g., Ref [1]).A lattice site and the edges around it are denoted by a star, as illustrated in Fig. 1a.The gauge operator corresponding to each star and group element g ∈ G is defined by: where s is the unitary operator corresponding to g applied to the gauge DoF on direction ê of s and the matter DoF on the lattice site s, respectively.In the nonabelian case the first pair of factors in Eq. ( 1) should correspond to a left group action, while the rest should correspond to a right action.We will focus on pure gauge models, and the matter charges will therefore be omitted from now on.As an example, we refer to gauge DoF with a Z 2 -symmetry, that is, the edges are occupied by 1/2-spins.The only nontrivial element in the group corresponds to the z Pauli matrix, σ z , and the gauge operator becomes In LGT models, the gauge operators on all of the stars in the lattice commute with the Hamiltonian and are therefore conserved.Gauge invariance requires that for any star s, in the absence of static charges, The requirement above imposes that the value of the gauge field going in each node (from the left and from below) must be equal to the value of the gauge field going out of the node (from the right and from above).For Z 2 gauge fields, this means that the number of downwardpointing spins around each lattice site must be even.
When G is a compact Lie group, the above may be intereperted as the Gauss law, by thinking of the nontrivial gauge charges as a flux: the flux going in each lattice site must be equal to the flux coming out of it, which enforces the fact that any excitation of the gauge sites can only exist on sets of edges that construct closed loops (see Fig. 1b), corresponding to the continuous Gauss law in the absence of charges, ⃗ ∇• ⃗ E = 0.The flux interpretation may also be valid in finite groups by introducing modular operators as in Eq. ( 2).In the Z 2 case, one may refer to the gauge upward-pointing spins as the vacuum, and downward-pointing spins as the flux.Downwardpointing spins may only appear in closed loops.
As carriers of interaction, it is relevant to discuss confinement and deconfinement in the context of models with gauge fields.In models with matter, the confinement-deconfinement phase transition relates to the strength of interaction between matter fields as a function of the distance between them, that is, the length of possible gauge-field-mediated paths between the matter particles.In pure gauge models, the above reduces into studying loops of gauge field excitations, also known as Wilson loops [2], i.e., operators that excite the gauge field around a closed loop of edges.In the Z 2 case, a Wilson loop is defined as where e runs over all of the edges that the loop W is composed of and σ x e is the x Pauli matrix acting on the spin on edge e.The expectation value of Wilson loops may decay as a function of their area (bulk size), in which case the state is in a confined phase, or as a function of their perimeter (boundary size), in which case the phase is deconfined.The confined phase is a disordered phase, and one may expect that it would imply less entanglement between parts of the system, since large flux loops have a small amplitude, while the deconfined phase is ordered and may exhibit larger entanglement.However, the relation of confinement and entanglement is not always straightforward, as we demonstrate and discuss below.Matter fields are placed on the lattice sites, and the gauge fields, mitigating the interaction of matter, are placed on the edges of the lattice.In yellow a star is denoted, that is, a matter site and the gauge sites it interacts with, and in blue a plaquette, that is, the minimal structure of a closed loop on the lattice.(b) In the pure gauge case, we omit the matter sites.In this case, gauge invariance is expressed by the requirement that the excitation on the gauge fields may only exist on edges that form closed loops.

A. Projected Entangled Pair States (PEPS)
In this work we focus on PEPS [93][94][95], a TN ansatz used in the representation of many-body quantum states.It is specifically relevant to states that exhibit area-law entanglement [96], that is, states in which the entanglement between two subsystems scales with the size of the boundary between them.
Consider a system of N lattice sites, with Hilbert space dimension d per site.A PEPS representation of the system's state will be composed of N tensors corresponding to the N lattice sites.Each tensor has one entry of dimension d, denoted as the physical leg, representing a Hilbert space that corresponds to the Hilbert space of a single site, and additional entries called the virtual legs, which connect tensors of nearest-neighbor sites.The di-mension of the Hilbert space represented by the virtual legs is denoted by D and referred to as the bond dimension.Contracting the virtual legs of all nearest-neighbor sites will result in a d ⊗N tensor representing the state of the system, which may be reshaped into a d N vector, the standard representation of a quantum state.The above is demonstrated in Fig. 2a.
A family of algorithms relying on this construction called infinite PEPS (iPEPS) is noteworthy: The system is assumed to be infinite and translationally invariant, and therefore defined by a repeating PEPS tensor.While contracting an iPEPS has been shown to be a computationally hard problem [92,97], an approximate contraction may be done, for example, using the boundary matrix product states method [98] or the corner transfer matrix method [99][100][101], allowing for the computation of expectation values and other properties of the state represented by an iPEPS.The PEPS and iPEPS ansätze are central numerical tools in the study of strongly correlated two-dimensional systems, and has been used for finding ground states [102][103][104][105][106][107][108], thermal states [109][110][111][112][113][114] and non-equilibrium steady-states [115][116][117][118][119][120][121] in two spatial dimensions.
When discussing density matrices rather than states, as is often the case in the study of entanglement, it is useful to define projected entangled pair operators (PEPO), and the operator version of iPEPS, iPEPO.A PEPO node has two physical indices rather than one -corresponding to two physical entries of the density matrix.A PEPS tensor may be used to create a PEPO tensor analogously to the definition of a density matrix out of (a) A state of an N -sized system with Hilbert dimension size d per site represented as a PEPS.The representation is composed of N tensors with a physical leg (pink) of dimension d and four virtual legs (blue) of dimension D. Contracting the virtual legs would result in a d ⊗N tensor, which may be reshaped into a d N vector.(b) The entanglement between two parts of the system is upper bounded by the total dimension of the virtual legs that cross the boundary between the two systems.(c) A PEPO tensor may be constructed by taking two copies of the PEPS tensor, performing a complex conjugation to one of them, and placing them back to back.(d) A single tensor in the gauge-invariant PEPS representation corresponds to a single lattice site and the edges on top of it and to its right, with two physical legs of dimension d in the pure gauge case.
a state vector, ρ = |ψ⟩ ⟨ψ|.The PEPO tensor is constructed by taking two copies of the tensor, performing a complex conjugation to one of them, and placing them back to back, as demonstrated in Fig. 2c.

Entanglement in Tensor Network States
We start by defining the entanglement measures that will be referred to in this work: For a given system in a pure state |ψ⟩, composed of two subsystems A, B, the reduced density matrix (RDM) of A is defined to be where Tr B denotes the tracing out of the degrees of freedom of subsystem B. The entanglement between subsystems A and B is standardly quantified by the von Neumann (vN) entropy [122]: Complementary entanglement measures are the Rényi entropies [122]: Note that lim n→1 S (n) (ρ A ) = S(ρ A ).While the Rényi entropies are entanglement monotones, they do not possess some useful properties of the vN entropy [122].However, the fact that their computation does not require taking a log of the matrix but only multiplying n copies of it, makes it more accessible numerically and experimentally [123][124][125][126][127][128][129][130][131][132][133][134][135], and sometimes also useful analytically by means of the replica trick, as can be seen, e.g., in Refs.[43,44,136] and in Sec.V A below.
We study the relation of the bond dimension D to the entanglement between two subsystems A, B of the system: All quantum states may be written in their Schmidt decomposition [137]: where {|i⟩ A } , {|i⟩ B } are orthonormal sets of states in the Hilbert space of A, B, respectively.From this decomposition, it straightforwardly follows that ρ A , ρ B have the same eigenvalues, . D is thus the number of nonzero eigenvalues (rank) of the RDMs.Any TN state may be brought to a form in which each set of bond indices residing on the boundary between subsystems A and B corresponds to a single coefficient ψ i [137], often referred to as the orthogonal or isometric form [138].In a PEPS, D is bounded from above by the dimension of all of the virtual legs crossing the boundary between A and B, D ≤ D |boundary| , as illustrated in Fig. 2b.Therefore, it is clear that the requirement that D is constant (as in translationally-invariant representations such as iPEPS) or efficiently dependent on the system size (as in all efficiently-representable PEPS states) implies that the state's entanglement must obey an area-law for any partition A, B.

B. Gauge-Invariant PEPS
In this section we briefly review the construction of gauge-invariant PEPS, and focus on Abelian gauge symmetries.For a detailed review, including the generalization to nonabelian symmetries, see Ref. [1].Based on the interpretation to the Gauge field excitations as flux, the PEPS tensor is chosen to represent a single lattice site and the edges coming out of it to the top and to the right (see Eq. ( 10)): The flux going into the node from the bottom and left must be equal to the flux going out of the node from the top and right.This may be imposed by assigning the Hilbert space corresponding to the virtual legs that has the same symmetry as the one on the physical gauge sites, as detailed below.
The local operators applied to each gauge link are elements of some finite or compact symmetry group G.The unitary operator corresponding to the group element g ∈ G is denoted by Θ (g) , and obeys where j corresponds to states in the Hilbert space of a single gauge link (group representations) and ϕ (g) j is the acquired phase.
We utilize the construction above and span the Hilbert space on the virtual legs using irreps of the same group G.The Hilbert space corresponding to the virtual legs may be of different dimension, but the symmetry of G must be obeyed.In this way, the flux value on the virtual legs becomes equivalent to the virtual flux value on the physical gauge legs.As in the physical Hilbert space, we define the unitary operator corresponding to the group element g ∈ G on the bond Hilbert space to be θ (g) : We now impose the local symmetry, as required in Eq. (3).The gauge operators apply to the physical legs of three neighboring PEPS tensors: The requirement in Eq. ( 8) may be translated into a connection between the physical and virtual Hilbert spaces: In order to impose Eq. ( 9), we require that the only nonzero terms in the gauge invariant PEPS tensors are the ones in which the outgoing flux on the top and right physical legs equals the outgoing flux on the top and right virtual legs, respectively.The total flux on the outgoing top and right virtual legs is required to be equal to the total flux on the ingoing bottom and left virtual legs: where .
= stands for equivalence of the physical and virtual group elements.
The states represented by PEPS that obey the requirement above are gauge-invariant by construction.Since the values of the physical indices are repeated on the upper and right virtual legs, we omit the drawing of the physical indices below for simplicity.Ref. [139] relies on the construction above to find PEPS representations to the ground state of the toric code model with added local magnetic field, and obtain an ansatz that requires D = 4 for the confined phase and D = 2 for the deconfined phase, as one might expect, recalling that the larger the bond dimension, the larger the entanglement may be.Ref. [140] describes a minimal model for Z 2 gauge DoF with rotational invariance with D = 2: (11) and analyze it analytically for several cases.Specifically, when β ≪ δ ≲ α, the PEPS is shown to represent a confined state if γ = 0, where γ > 0 results in a deconfined state.
One may also construct PEPS-representable states with an infinite symmetry gauge group, such as U(1), by truncating the link Hilbert space in the representation basis, and keeping only a finite number of representations [141] (for example, in U(1) where representations corresponds to the angular momentum value along the group circle, one may restrict the link Hilbert space to angular momenta −m . . .m, resulting in d = 2m + 1).The physical Hilbert space dimension d then becomes finite, and therefore the bond dimension may also be finite and the state representable by PEPS.

IV. SYMMETRY-RESOLVED ENTANGLEMENT IN GAUGE-INVARIANT STATES
Before we dive into the many definitions in this section, it is noteworthy that in gauge-invariant theories, the structure of the Hilbert space is different from the general case: In a general system, the Hilbert space of systems A, B is of the structure However, once gauge-noninvariant states are excluded from H AB , the Hilbert space loses its tensor-product form, and the Hilbert spaces H A , H B become ill-defined, and in turn also ρ A , ρ B and the entanglement measures which rely on them.In this paper, we follow the standard method for solving this problem [34,[142][143][144][145] by embedding the gauge-invariant Hilbert space in the full Hilbert space H A ⊗H B .As we will shortly explain, when defining symmetry-resolved entanglement and limiting the Hilbert space into a single block in ρ A and its corresponding block in ρ B , the problem above fades away: The gauge-invariant Hilbert subspace is naturally in a product form, and the embedding becomes trivial.We will continue concentrating on the Abelian case, and will briefly consider the nonabelian case in Sec.IV B.
First, we show that thanks to the Gauss law, the RDM ρ A is block diagonal, with different blocks corresponding to different values of the flux going through the boundary between subsystems A, B, as illustrated in Fig. 3b.We work in the eigenbasis of the unitaries Θ g , {|j⟩} d j=1 , and write the state of the full system as where ⃗ j e is the state on the site e.For each star s on the boundary of subsystems A, B, we denote by s A , s B the star parts residing on subsystems A, B, so that s = s A s B .
From Eq. ( 3) it is required that the value of the total flux on s A fully determines the value of the total flux on s B .We denote by ϕ s , ϕ s the total flux values on the star parts s A , s B across the boundary in subsystems A, B, respectively, and by ⃗ j ϕs , ⃗ j ϕs the subspaces of A, B that correspond to the flux ϕ s .The state may then be written as Tracing out the DoF of subsystem B, and noting that states corresponding to different fluxes are orthogonal, the remaining RDM is The obtained RDM is block diagonal, each block corresponding to a different value of ϕ s .One may repeat this argument for all stars on the boundary, and obtain a block diagonal RDM where each block corresponds to the set of flux values {ϕ s } s∈boundary[AB] , which for brevity we will henceforth denote by ϕ.We denote by ρ A the block corresponding to the flux combination ϕ, and note that ρ A = ⊕ ϕ ρ ϕ A .One may now define the symmetry-resolved entanglement as the entanglement obtained from each block separately: the symmetry-resolved vN and Rényi entropies are defined in analogy with Eqs.(5,6), respectively:

S
(n) Note that the symmetry-resolved vN and the moments of the Rényi entropies sum to their nonresolved analogs, . It is also interesting to define the normalized symmetry-resolved entanglement measures, obtained by normalizing the blocks,

S
(n) We follow common notation and denote the different blocks ρ A (ϕ) by superselectionresolved (SR) entanglement.In the rest of this paper, we follow the normalized definition in Eqs.(14,15), as justified later in this section.The case in which the symmetry-resolved entanglement is independent of the block is called equipartition [45].It is observed analytically and numerically in many systems with one and two space dimensions, although some models without equipartition have been observed in 1+1D systems [68,76,[146][147][148][149][150].We now note an important observation: The contribution to the full (nonresolved) entanglement comes from two sources: the division of ρ A into blocks and the contribution of each block (SR entanglement).For example, for the von Neumann entropy one has where p(ϕ) = Tr ρ (ϕ) A . The off-diagonal terms (coherence) in ρ A is limited to the superselection sector blocks, that is, relevant to the first term only.If one limits themselves only to gauge-conserving operators and observables on the stars and star-parts in subsystem A (whether equipartition applies or not), the contribution of ρ A 's division into blocks becomes inaccessible: The effect of the division into blocks on any observed quantity can be simulated classically.This makes the normalized SR entanglement a measure for the accessible entanglement, sometimes referred to as distillable entanglement.The limitation to gauge-conserving operators applies in simulated high-energy models, and may also be applied in condensed matter models, further motivating the study of SR entanglement in such gauge-invariant states.

A. Superselection-resolved entanglement in gauge invariant tensor network states
In this section we derive our main result, characterizing the behavior of superselection-resolved entanglement in gauge-invariant lattice states.As in Sec.III, we focus on a pure gauge model on a square lattice obeying an Abelian symmetry; the nonabelian case will considered in Sec.IV B. The generalization to different lattice types, higher dimensions, or gauge models with matter is straightforward.
Our analysis is based on the relation of entanglement and bond dimension described in Sec.III A. This relation has two implications: (a) The projection to a specific superselection sector, that is, a specific set of RDM eigenvalues, may be obtained solely by applying a projection to the virtual Hilbert subspace, as illustrated in Fig. 4.Such a projection becomes natural for PEPS constructed as in the protocol described in Sec.III B: A projection onto a specific physical charge turns into a projection onto the corresponding charge on the virtual legs.If a star-part contains only one of the lattice sites in a PEPS tensor, the decomposition of the tensor is required, as depicted in Fig. 4b.(b) After the projection onto a superselection sector on the boundary, the dimension product of the virtual indices bounds the number of nonzero eigenvalues of ρ A (ϕ) , ρ B (ϕ) , which in turn bound the SR entanglement from above.This understanding is the basis of our results.
The projection onto a single block is obtained by applying many local projections, on all of the star-parts on the system, see Fig. 4d.Note that the partition boundary divides the stars into subgroups of three and one gauge sites (Fig 4a,b).Determining the flux on a star-part on the edge therefore fully determines the flux on the single-site star-part.On the corners of the partition, the star-parts on both subsystems contain two gauge sites.Determining the flux on one of the star-parts therefore does not determine the flux of any single gauge site, and therefore also does not determine the charge on any single virtual leg (Fig. 4c).
We refer to the case where D = d.This requirement is obeyed in PEPS representations of various models, such as the toric code model with string tension, displaying a topological-polarized phase transition [151], thermal state inspired TN states undergoing a symmetry enriched topological phase transition [152], and the TN construction proposed in Ref. [140] and presented in Eq. ( 11), displaying both confined and deconfined phases.We recall that making a projection of a single site into a specific charge is equivalent to projecting the virtual leg next to it onto the corresponding charge.In the case of D = d, this results in a projection of the virtual leg on the boundary between the subsystems onto a Hilbert space of size 1.Effectively, the projected tensor has bond dimension 1 on all of the edges on the subsystems' boundary, adding no contribution to the bond dimension product and in turn, to the SR entanglement.The only exception are the corners of the system, which may require a bond dimension D = d on the corner.We then obtain a dependence of the entanglement on the number of corners on the boundary between subsystems A and B -a corner-law entanglement: which to our knowledge, has not been observed in any other system before.We note that one may choose the bipartition around the corners such that the star part sizes are odd and the corner contribution vanishes, as can be seen in Sec.V B below, but this bipartition may be considered less natural as the one in Fig. 4d.
We stress again the point made in Sec.IV: If one restricts themselves to applying and measuring only gaugeinvariant operators, as is often the case, the only quantum correlation between the systems is the SR entanglement.The correlations between two parts of the system stemming from the division into superselection sectors is completely classical, in the sense that it can be simulated classically and may not be harvested, for example, to create Bell pairs.Our result shows that for all states in which D = d, the SR entanglement is constant and small for typical partitions (e.g., rectangular systems), and the quantum correlation between the subsystems in this case is thus strictly limited.
In the general case where D > d, the projection in Fig. 4 results in virtual legs of dimension D/d.The SR entanglement then obeys an area law, proportionate to the number of virtual legs on the boundary between subsystems A and B: One may recall again the result of Emonts et al. [139]: A deconfined phase, intuitively believed to have larger entanglement, typically requires D = 2d, where the confined phase may be represented to a good approximation using D = d.The relation between entanglement, whether SR or not, and confinement, is not always intuitive as is the case in Ref. [139].For example, the model in Eq. ( 11) is constructed to have D = d, but may represent both confined and deconfined phases.In appendix A we present the "perfectly confined" state in d = 2, i.e., a state for which a Wilson loops' expectation value on its area is exactly exponentially-dependent on its area.Such a state requires D = 4, but as we show in appendix A, its SR entanglement is constantthe entanglement's dependence on system's properties is trivial and not even corner-law.Ref. [153] shows that in an isometric mapping between Kitaev's toric code model and a LGT model, both presenting a confined-deconfined transition, the entanglement spectrum changes, indicating that entanglement in gauge-invariant models indeed behaves differently in LGT than in general models.From the above, we understand that the relation of entanglement and confinement in LGT models requires further study.The physical charge on the star part is denoted by q.Since the total charge on the star must equal 0, the charge on the bond must equal the bond-representation equivalent of −q.

𝑃
(b) When the completing site of the star belongs to the same PEPS node as one of the sites on the star, the node may be split into two tensors, and the projection applied to the new bond dimension.(c) In the case of a two-site star part (corner), the projection is done on two virtual legs together and imposes the total charge on both of them, rather than the chrage of each of them separately.(d) The projection onto a single superselection block is done by projecting all star parts.

B. Generalization to nonabelian theories
In nonabelian theories, the requirement (8) becomes slightly more complex and therefore also the PEPS representation of gauge-invariant states.However, the idea behind our argument remains the same, and a corner law is obtained in the representation with minimal bond dimension, as we now explain.
In nonabelian symmetries, the gauge requirement Eq. ( 3) becomes a requirement that all stars are singlets of the gauge symmetry group.The states of each physical leg are then of the form |jmm ′ ⟩ , where j stands for the representation of the group and m, m ′ are states in the representation j, corresponding to left-and right-applied operators (for full details, see Ref. [1]).The bond states are then of the form j m , where j is a group representation and m an element in j.
To form a singlet, the two parts of a star should correspond to some representation j and its conjugate, and the values of m should match as well.The density matrix of A then breaks into blocks corresponding to specifying j and m for all the stars cut by the boundary of A. As a result, when the bond dimension is minimal, so that each representation in the physical space appears once in the virtual space (in which case the virtual bond dimension is smaller than the physical one, since the virtual bonds carry a single index m, while the physical bonds carry two indices m, m ′ ), for each block all the virtual legs at the boundary are fixed except for the corners, leading again to a corner law.
We also note that, for each representation j and its conjugate the amplitude of the different states m in the singlet state is equal in magnitude, leading to equipartition between blocks with the same j value [39,145]; equipartition between blocks with different j values is not guaranteed by symmetry, as in the Abelian case.Finally, As in the Abelian case, gauge models with infinite symmetry groups such as SU(2) may be truncated keeping a finite number of representations per link, such that they are representable by PEPS and the argument applies to them as well.

V. NUMERICAL RESULTS
A. Transfer matrix method for system-size dependent properties In this section, we go over the technique we used for computing the confinement measure, i.e., area-(bulk-) law dependence of Wilson loop expectation values, following Ref.[140].Using the same line of thought, we follow Ref. [136] and describe the method we used for computing the area-(boundary-) law dependence of full and SR Rényi entropies.
Wilson loop expectation values may be computed by a tiling of PEPS tensors.We denote traced PEPO nodes, i.e., PEPO nodes with contracted physical indices, as follows: The expectation value of an R 1 × R 2 Wilson loop by the multiplication of matrices constructed of rows of traced PEPO tensors: Based on Eq. ( 20), we define the transfer matrices composed of contracted rows of PEPO tensors: We briefly analyze the dependence of Ŵ on R 1 , R 2 as is done in Ref. [140].The largest eigenvalues of R1 due to the repetition of E || and E in the numerator and denominator, respectively.By requiring consistency of the computation when rotating the system by π/2, Ref. [140] obtains an identical dependence of Ŵ on R 1 .
Ref. [136] shows that the dependence of r ′ 1 (R) must be exponential: where Γ, κ are constants.κ > 0 would result in an arealaw contribution to Ŵ -confinement.Therefore, we define κ to be the measure of confinement in our numerical computations below.r ′ 1 (R) is computed using the TEBD method for non-unitary evolutions [154][155][156] for several values of R, and κ is obtained from fitting the data as a function of R.
Ref. [136] brings up a similar mechanism for the study of Rényi entropies, and specifically, the second RDM moment also known as the purity.The computation of the purity of an R 1 × R 2 rectangular system may be done by using two copies of the PEPS density matrix.A single row of the TN construction, analogously to the Wilson loop computation, would be: (24) One may compute the SR-resolved purity of the block containing no flux on all star-parts.If the reasonable assumption of entanglement equipartition applies, p 2 ({↑, ↑ . . .}) is equal to all other SR purities.The SR purities are obtained by applying the local projection on the bond indices: Following Ref. [140], we perform a similar analysis as in the Wilson loop case above.The largest eigenvalues of Ẽ, Ẽ|| (R 2 ) are denoted by r1 , r′ 1 (R 2 ).In the entanglement case, the PEPS construction imposes that p 2 , p 2 ({↑, ↑, . . .}) has only area (boundary)-law contributions, therefore it is deduced that r′ 1 (R) = const.always.We quantify the area-law by (r ′ 1 /r 1 ), and expect that in the case of a corner-law entanglement, (r ′ 1 /r 1 ) = 1.We denote as the area-law measure.

B. Numerical results
We start by demonstrating the corner-law entanglement by studying a D > d model.We simulate a pure Z 2 gauge model with D = 4.The iPEPS nodes obey the restriction in Eq. (10), and the nonzero elements are chosen independently and randomly from a Gaussian distribution N (µ, σ).We fix the mean µ = 1.When σ = 0, Kitaev's toric code ground state is obtained: where p stands for plaquettes of the lattice as demonstrated in Fig. 1, X p = ⊗ e∈p σ x e , and N is the size of the system.This model may be represented by D = d = 2 by choosing α = β = γ = δ = 1 in the model defined in Eq. (11), and is therefore expected to have no area-law contribution to the SR entanglement; only corner-law contribution is allowed [157].As σ is increased, the structure is expected to be lost, and an area-law SR entanglement is observed.We compute the dependence of the SR and full purities on system size as explained in Sec.V A. In Fig. 6a, η 2 (R) is presented for the full and SR purity for an open system of width 2R where R = 14, as computed by averaging over 10 occurrences of the random iPEPS tensors.Indeed, the area-law measure η 2 of Eq. ( 26) vanishes in the SR purity when σ = 0, and becomes larger as σ increases.As a result, the full entanglement's dependence on area also grows, where when σ = 0, the expected value η 2 = 1 is obtained.
We then move to study a model with d = D, as defined in Eq. (11).We choose the parameters β ≪ δ ≲ α, in which, following Ref.[140], we expect to see a confined phase at γ = 0 and deconfined phase at γ > 0. Specifically, we choose α = 1, β = 0.3, δ = 0.9 and 0 ≤ γ ≤ 2. First, the confinement of the phase is computed as explained in Sec.V A. r ′ 1 (R) of Eq. ( 22) is computed for E || (R) (Eq.( 21)), as defined in an open system of 64 sites, for R = 14 . . .32. κ is then extracted and plotted in Fig. 6b.As expected, when γ ≪ 1, κ > 0 and the phase is confined, and κ = 0 deconfinement is observed for all other values of γ.As computed in Ref. [140], we obtain κ As a sanity check, we compute the dependence of p 2, p 2 ({↑, ↑ . . .}) of rectangular systems for R = 6 . . .14, and compute the area-law SR entanglement as discussed in Sec.V A. The results are plotted in Fig. 6b.As expected, there is no dependence of the SR entanglement on the area in any value of γ.
We then study the SR entanglement as a function of the number of corners in the system.We study a stairslike system geometry which allows for a constant area (boundary) length for a different number of corners, as illustrated in Fig. 5.Note that the bipartition is chosen such that not all corners contribute to the entanglement, as can beseen in Fig. 5.In our computation, the height and width of the system is L = 6, and the number of corners is c = 1 . . .6 accordingly.The computation was performed with the BMPS method [98], simulating an infinite system with open boundary conditions.The SR purity is computed for two blocks -The block corresponding to s A = 1 in all star-parts, ϕ = {↑, ↑ . . .}, and one block with a random flux distribution on the star parts.Since d = D in this model, we expect that where b 0 , b 1 are constants.In Fig. 6b, we plot the cornerlaw measure b 1 for both blocks.Both blocks display a virtually-identical behavior, supporting the applicability of equipartition in this model.It is interesting to note that the behavior of the entanglement, displaying a large entanglement around γ = 1 and small entanglement (SR and full) far from it.The dependence of the SR entanglement on γ is in no correlation with confinement in this model.As mentioned above, the relation of SR entanglement to confinement (or lack thereof) requires further study.

𝐿 𝐿
. The geometry we used for studying the cornerdependence of a D = d model.The number of corners (three in the illustration) may be changed while keeping the area (boundary) of the system constant.Note that the bipartition was chosen such that all corners but the bottom-left ones divide the stars into odd-sized star parts, and therefore they do not contribute to the entanglement.

VI. CONCLUSION AND FUTURE OUTLOOK
We have used the tensor network formalism to study superselection-resoved (symmetry-resolved) entanglement in gauge invariant models.SR entanglement is particularly interesting in such models, since it is the only accessible (distillable) entanglement when the gauge-invariance is strictly imposed in a system's partition.We have shown that when the state is representable by PEPS with D = d, which is the case in a variation of interesting states, the SR entanglement is bounded by a corner-law, which implies a constant entanglement in  11).The confinement measure κ is presented in the top panel, the full entanglement's area-law dependence is in the middle panel.and in the bottom panel we present the area-and corner-law dependence of the SR entanglement.Note that the behavior of the SR entanglement is very similar, consistent with equipartition.
most interesting geometries (such as rectangular or cylindrical subsystems).
Our argument may be readily used for different lattice types, nonabelian gauge symmetries, higher dimensions (in which an edge-law is obtained), and models with matter.However, the variety of effects and phases in the cases mentioned above has not been thoroughly studied.It would be interesting to explore the range of phases and models representable by D = d and the implications of our analysis on such systems.
Lastly, as we see in several past studies as well as our own, it seems that the relation between confinement and entanglement, SR or general, requires further examination.We hope that our work contributes to future studies of this relation.
Our  In this section, we discuss the relation between the SR entanglement and confinement by examining the SR entanglement in states that display an exact area-or perimeter-law decay of Wilson loop expectation value.We show that in both cases, the SR entanglement converges to a constant value in the limit of large subsystems.From this, it is apparent that non-trivial SR entanglement comes from the behavior of the Wilson loop expectation values in finite loop size, or from effects unrelated to flux and confinement.It would be interesting to fully characterize the relation of SR entanglement and confinement.The directionality of the loop in the environment determines which plaquettes in the system will contribute to the amplitude: In yellow, the loop is closed from below the system and a single plaquette contributes, and in pink, the loop is closed from above and three plaquettes contribute.

X X X X
We start with a perfectly-confined phase, i.e., states in which the Wilson loops expectation value is exactly exponentially-dependent in the Wilson loop area: where Γ a , κ a are constants.We observe a perfectlyconfined state of a Z 2 pure-gauge field system: where p counts the plaquettes in the system as demonstrated in Fig. 1 in the main text, and X p = ⊗ e∈p σ x e .For simplicity, we study the SR entanglement in blocks with a single flux line going in and out of them, as in Fig. 7, but the generalization to all blocks is straightforward.The state in Eq. (A1), projected into the block denoted by {q}, may be written in the computational basis as where A(i, j) stands for the total area of the flux loops in the total system in the state |i⟩ B |j⟩ A .From Fig. 7, we conclude that the area of the loop is determined based on the directionality of the closing of the loop: If the loop is closed below the subsystem (yellow route in Fig. 7), the loop area would be the sum of the area covered by the loop in the environment and the area below the flux line in subsystem A. If the loop is closed above the subsystem (pink route in Fig. 7), the total area would be the sum of the loop area in the environment and the area above the flux line in A.
The projected state of the system may thus be written as where A(A) is the area of subsystem A, that is, the number of plaquettes in A. Note that adding inner loops in A or in B is consistent with the above description.
We now denote The obtained RDM is enough to show that the SR entanglement is independent on subsystem A's size and features: Since ρ A (ϕ) may be written as a 2 × 2 matrix, the maximal number of nonzero eigenvalues is 2, bounding the entanglement from above by [S A (ϕ)] max = log 2. However, for completeness, we continue the analysis of ρ A (ϕ) below.
In the basis {|0⟩ A , |1⟩ A } , the block ρ A (ϕ) may be written as where are the contributions of the environment to the states' amplitude for flux lines closed above and below subsystem A. Assuming an infinite system, one may assume C b = C t (this assumption is exact for a system a reflection symmetry).The block ρ A (ϕ) then becomes The SR vN entropy would then be , where the approximation is applicable due to the exponential decay of v 0 with subsystem A's size, as can be seen in Eq. (A2).The entanglement therefore converges to a constant (and small) entanglement as A's size becomes larger.Note that the above analysis is identical whether the boundary conditions of the system are open or periodic.
Secondly, we study a perfectly-deconfined phase, i.e., states with an exact exponential-dependence of the Wilson loop expectation value in its perimeter: where Γ p , κ p are constants.Such a state may be obtained by adding string tension to Kitaev's toric code model, as in Refs.[151,[158][159][160], where e runs over all of the edges (gauge sites) in the system.Here, one may quickly see that the SR entanglement between two subsystems vanishes: The state of the system, projected into a specific block {q}, may be written as which is a product state, disentangled by definition.We note in passing that our results could be easily extended to a model wave function with nontrivial values of both κ a and κ p , displaying a confined phase and small SR entanglement as in Eq. (A1).

Figure 1 .
Figure1.The discretization of a gauge model on a lattice graph.(a) Matter fields are placed on the lattice sites, and the gauge fields, mitigating the interaction of matter, are placed on the edges of the lattice.In yellow a star is denoted, that is, a matter site and the gauge sites it interacts with, and in blue a plaquette, that is, the minimal structure of a closed loop on the lattice.(b) In the pure gauge case, we omit the matter sites.In this case, gauge invariance is expressed by the requirement that the excitation on the gauge fields may only exist on edges that form closed loops.
Figure 2.(a) A state of an N -sized system with Hilbert dimension size d per site represented as a PEPS.The representation is composed of N tensors with a physical leg (pink) of dimension d and four virtual legs (blue) of dimension D. Contracting the virtual legs would result in a d ⊗N tensor, which may be reshaped into a d N vector.(b) The entanglement between two parts of the system is upper bounded by the total dimension of the virtual legs that cross the boundary between the two systems.(c) A PEPO tensor may be constructed by taking two copies of the PEPS tensor, performing a complex conjugation to one of them, and placing them back to back.(d) A single tensor in the gauge-invariant PEPS representation corresponds to a single lattice site and the edges on top of it and to its right, with two physical legs of dimension d in the pure gauge case.
(ϕ) A by superselection sectors, and the entanglement measures S A (ϕ) , S

Figure 3 .
Figure 3. (a) The superselection sectors are subspaces of the Hilbert space of subsystem A, defined by the flux on the starparts on the boundary between the subsystem and its environment.In circles are examples of star parts on the boundary between the subsystems.(b) The RDM is deconposed of blocks, each block corresponding to a single superselection sector.

Figure 4 .
Figure 4.The projection onto specific flux of a star part may be done by projecting the virtual legs going into or out of the star part.(a) When the completing site of the star belongs to a different PEPS site than the sites in the star part, the projection may be done on one of the ingoing virtual legs.The physical charge on the star part is denoted by q.Since the total charge on the star must equal 0, the charge on the bond must equal the bond-representation equivalent of −q.(b) When the completing site of the star belongs to the same PEPS node as one of the sites on the star, the node may be split into two tensors, and the projection applied to the new bond dimension.(c) In the case of a two-site star part (corner), the projection is done on two virtual legs together and imposes the total charge on both of them, rather than the chrage of each of them separately.(d) The projection onto a single superselection block is done by projecting all star parts.

Figure 6 .
Figure 6.Numerical study of full and SR entanglement in pure Z2 gauge models.(a) Random tensors with D = 2d = 4.The SR entanglement's dependence on the area (boundary) vanishes when σ = 0 and the state may be represented with D = d, and grows as the tensor becomes far from it.(b) The generic D = d = 2 model, Eq. (11).The confinement measure κ is presented in the top panel, the full entanglement's area-law dependence is in the middle panel.and in the bottom panel we present the area-and corner-law dependence of the SR entanglement.Note that the behavior of the SR entanglement is very similar, consistent with equipartition.
Appendix A: Superselection-resolved entanglement in perfectly confined and deconfined states

Figure 7 .
Figure 7.The studied subsystem has two ingoing flux lines.The directionality of the loop in the environment determines which plaquettes in the system will contribute to the amplitude: In yellow, the loop is closed from below the system and a single plaquette contributes, and in pink, the loop is closed from above and three plaquettes contribute.
work has been supported by the Israel Science Foundation (ISF) and the Directorate for Defense Research and Development (DDR&D) grant No. 3427/21 and by the US-Israel Binational Science Foundation (BSF) Grant No. 2020072.NF is supported by the Azrieli Foundation Fellows program.EZ acknowledges the support of the Israel Science Foundation (grant No. 523/20).JK is supported by the Israel Academy of Sciences and Humanities & Council for Higher Education Excellence Fellowship Program for International Postdoctoral Researchers.