Entanglement and confinement in lattice gauge theory tensor networks

We develop a transfer operator approach for the calculation of R\'enyi entanglement entropies in arbitrary (i.e. Abelian and non-Abelian) pure lattice gauge theory projected entangled pair states in 2+1 dimensions. It is explicitly shown how the long-range behavior of these quantities gives rise to an entanglement area law in both the thermodynamic limit and in the continuum. We numerically demonstrate the applicability of our method to the $Z_2$ lattice gauge theory and relate some entanglement properties to the confinement-deconfinement transition therein. We provide evidence that R\'enyi entanglement entropies in certain cases do not provide a complete probe of (de)confinement properties compared to Wilson loop expectation values as other genuine (nonlocal) observables.


Introduction and motivation
Characterization of entanglement properties in quantum field and gauge theories is of upmost importance for the understanding of physical systems in condensed matter, particle and gravitational physics [1][2][3].As such, it allows to gain insights into important properties of emergent phenomena.In this article, we are particularly interested in the interplay of entanglement and confinement.The latter is a nonperturbative phenomenon, which describes the binding, i.e. confinement, of static charges.It appears prominently in quantum chromodynamics (QCD) as the theory of strong interactions in the standard model of particle physics, for which it manifests itself through the existence of color-neutral and strongly interacting quark bound states (e.g.hadrons and mesons) at low energies [4].On the other hand, at high energies, a deconfined quark-gluon plasma exists as a result of asymptotic freedom, for which a perturbative treatment can be possible [5][6][7].Examples of confinement occur also in lower-dimensional quantum many-body (QMB) systems and quantum field theories (QFTs) [8][9][10][11].They provide a way to study many long-standing open questions in this context using novel tools from a quantum information perspective.
One main motivation for our work comes from holography and the AdS/CFT correspondence [12][13][14].The seminal work [15,16] provided a holographic prescription for entanglement entropy, as the most prominent and rigorous entanglement measure of pure states, by showing that it can be calculated as the minimal surface area in the bulk belonging to a chosen boundary region.When studying this quantity, the authors of [17] proposed that entanglement (entropy) can serve as a probe of the confinement-deconfinement transition in gravitational duals of large-N c gauge theories.This idea was further studied and confirmed in many holographic models for QCD and beyond, see e.g.[18][19][20][21][22][23][24][25][26][27].Following this line of research, several recent works unveiled the importance of entanglement properties also directly in processes relevant for particle and nuclear physics [28][29][30][31][32][33][34][35][36][37][38][39].
In this work, we are using tensor network methods to analyze entanglement properties in (2+1)-dimensional lattice gauge theories (LGTs). 1 Tensor networks and tensor network states are a concept and algorithmic tool originating in quantum information science.They are formulated in a Hamiltonian framework and provide efficient ansätze for wavefunctions of QMB states with a polynomial number of parameters.(For broad reviews on that topic see [40][41][42][43][44][45][46].)In contrast to Monte-Carlo simulations in Euclidean spacetime, tensor networks are sign-problem-free and usable for real-time simulations.They also allow to implement gauge symmetries as the necessary basis for the study of gauge theories in particle and condensed matter physics [47,48].In one spatial dimension, matrix product states (MPS) have been successfully used to study a plethora of physical effects in spin chain models, QFTs and LGTs, see e.g. the recent reviews [44][45][46][47][48] and references therein.Specifically, they were also employed to study the effect of confinement on static and dynamical entanglement properties [11,[49][50][51][52][53][54][55][56][57][58].In two dimensions, projected entangled pair states (PEPS) are the natural generalization of MPS [40,42,44].Several gauging mechanisms have been proposed for them [59,60], which allow to study LGTs in (2+1)D.For example, the works [61,62] analyzed infinite systems using a translationally invariant (uniform) PEPS ansatz.The authors of [63][64][65][66] combined the gauging principle with a Gaussian fermionic PEPS, which can be efficiently contracted using Monte-Carlo methods [67].Beyond PEPS, also other Hamiltonian tensor network types [68][69][70][71][72][73] and non-Hamiltonian partition function approaches with tensor networks [74] have been employed for the study of higher-dimensional LGTs.
When studying entanglement properties for a bipartition of a physical system into a subsystem A and its complement B via entropic quantities, the reduced density matrix ρ A ≡ Tr B ρ is calculated by taking a partial trace of the full density matrix ρ over the Hilbert space of the complement region.This procedure assumes the separability of the underlying full Hilbert space into a direct product form, i.e.H A ⊗ H B .In LGTs, entanglement considerations are complicated by the fact that the space of physical states does not admit such a direct product structure: As a consequence of the Gauss law, which follows from demanding local symmetry constraints, there exist nonlocal degrees of freedom hindering this decomposition.Several approaches have been developed to overcome this obstacle and analyze entanglement properties in LGTs [75][76][77][78][79][80][81][82][83][84][85][86][87][88][89][90][91]. 2 In the present article we work within an extended Hilbert space approach [77,78,[81][82][83], which is obtained by taking the tensor product of all gauge field Hilbert spaces on the links of the lattice.A partial trace over the complement links is then well-defined.Alternatively, one can consider the entanglement for superselection sectors [79,80], which is underlying a symmetry-resolved entanglement approach [92,93] that is developed independently in the companion paper [94] for LGTs.
In pure LGTs, Wilson loops -operators creating closed flux lines -are important nonlocal, gauge-invariant observables.The decay of their expectation value allows to probe confinement properties: While a Wilson loop area law is present in the confined phase, a perimeter law scaling implies the existence of deconfined static charges.In the recent work [95] it was shown how the gauged PEPS setup introduced above can be employed to study Wilson loops and hence (de)confinement properties for arbitrary (i.e.Abelian and non-Abelian) pure LGTs.This construction was based on a transfer operator approach.Specifically, it was deduced how local properties of transfer operators, which are constructed out of the tensors obeying the gauge symmetry, dictate the long-range behavior of the Wilson loop.Here, we are following a parallel route to show how (normalized) Rényi entanglement entropies can be constructed for the two-dimensional LGT PEPS using transfer operators.We explicitly show how an entanglement area law [96] emerges in the thermodynamic limit from our general result and persists in the continuum. 3Based on these general results, we will argue that Rényi entanglement entropies are not an equally rigorous measure of confinement or deconfinement as the Wilson loop.However, using numerical techniques, we analyze how confinement properties do leave imprints on entanglement measures for the Z 2 LGT as a specific example.
We finally would like to add that our explorations are also motivated by recent advances of quantum simulation experiments [97][98][99], in which confinement properties such as meson physics have already been explored [100][101][102][103] or are at reach of current technologies [104].At the same time, quantum simulations also allow to measure and explore entanglement properties of quantum many-body systems and LGTs [105][106][107][108][109][110].In this line of research, we are addressing in this paper the timely topic of the interplay of both of these phenomena for two-dimensional LGTs, and want to provide further stimulus for its experimental exploration.
This paper is organized as follows.In section 2 we review the underlying gauged PEPS ansatz and elaborate on the construction of transfer operators.In section 3 we develop the transfer operator approach for the calculation of Rényi entanglement entropies.This idea is at first described in 1D for MPS and then generalized to 2D for PEPS.We subsequently calculate Rényi entanglement entropies by contracting transfer operator plaquettes within a 2D lattice and derive their long-range properties.In section 4 this formalism is numerically applied to the Z 2 LGT and contrasted with confinement properties therein.Concluding discussions and an outlook are given in section 5.The appendix A contains mathematical derivations of properties of a dominant transfer matrix eigenvalue.

Tensor network construction
In this section we review important concepts of LGTs and the gauged PEPS ansatz.We then outline the construction of transfer operators in this framework.We intend to focus on the most crucial properties, which are relevant for the following sections in this paper.For further mathematical background and a comprehensive discussion of LGTs themselves, we refer to [95], whose notation we will mainly follow in this section.

LGTs and gauge invariant PEPS
LGTs describe the interaction of matter particles through the mediation by gauge fields.In this paper we are studying (2+1)D LGTs in a Hamiltonian formalism [111], for which the spatial coordinates are discretized, while time is kept continuous.We restrict ourselves to the pure gauge case, in which the gauge fields are placed on the links of a 2D lattice, while fermionic matter fields, which would reside on the lattice sites, are absent; cf.Fig. 1.
Let us denote the lattice sites as x ∈ Z 2 , the directional unit vectors as êi , and the emanating links as l ≡ (x, i), where i = 1, 2 labels the horizontal and vertical direction, respectively.By definition, LGTs are invariant under a local, i.e. spacetime-dependent symmetry.We introduce a gauge group G and denote its group elements as g ∈ G.Each link l hosts a gauge field Hilbert space H l , which can be either spanned by the gauge group element states |g⟩, or, conventionally, by dual representation basis states |jmn⟩.For these, j labels an irreducible representation (irrep) of G, and m, n label identifiers within each.
Let us consider the gauge transformation For a compact Lie group G, for which the operators can be written as Θ g = exp(iϕ a (g)R a ) and Θg = exp(iϕ a (g)L a ) in terms of their right (R a ) and left transformation generators (L a ) and group element dependent parameters ϕ a (g), the gauge transformation (2.1) can be reexpressed using the Gauss law operator The gauge invariance condition (2.2) is then reformulated in terms of the Gauss laws Here, and in the following, we assume the absence of static charges.Since G a (x) can be seen as the divergence of electric field values, eq.(2.5) tells us that the ingoing flux at each site (from the left and down links) equals the outgoing one (on the right and upper links).
In this paper we make use of the gauged PEPS ansatz of [60], which allows us to study a pure LGT using a tensor network construction.(For further explorations of this formalism see [63,64,66,67].)We assume a 2D lattice with periodic boundary conditions in both directions.On each lattice site a tensor A st ruld is placed.Fig. 2(a) shows a graphical representation of this tensor using diagrammatic tensor network notation. 4The tensor A st ruld has 4 virtual indices (or legs), labelled as r, u, l, d for the respective directions right, up, left and down.The l and d legs are considered ingoing, while the r and u legs are the outgoing ones.They are associated with the virtual Hilbert spaces H r , H u , H l and H d , which are spanned by group multiplet states |jm⟩ (j labels irreps of G, m an identifier within each).The numerical size of these virtual indices is denoted as the bond dimension χ and allows to vary the number of variational parameters.While all states within a multiplet must be included, the overall number of multiplets may be truncated (for detailed discussions see [60,95]).The additional indices s and t (standing for side and top) correspond to the physical gauge field Hilbert spaces H s and H t , placed respectively on the r and u links.As discussed above, they are spanned by the dual representation basis states |jmn⟩, in which the irreps might be truncated as well, that is, a subset of the irreps j of the group is included.The resulting physical-virtual state |A(x)⟩ = A jsmsns;jtmtnt jrmr;jumu;j l m l ;j d m d |j s m s n s ; j t m t n t ⟩ |j r m r ; j u m u ; j l m l ; on each lattice site is therefore an element in the Hilbert space The second building block of the PEPS are projector states B 1,2 on the links (cf.Fig. 2(b)), defined as (2.7) Their purpose is to project the virtual states between two neighboring sites onto maximally entangled states via contraction.In this way, the contraction of all virtual tensor indices over the whole lattice, |Ψ⟩ = x,i yields a PEPS |Ψ⟩ with only physical (gauge) degrees of freedom, cf.Fig. 2(c).
We assume the tensor network to be translationally invariant, i.e. the same tensors are placed on all sites.In order for the resulting PEPS |Ψ⟩ to be gauge invariant, it has to obey the condition (2.2).Following the construction in [60], one can directly see that this condition is fulfilled (cf.Fig. 1) if the physical-virtual state satisfies (∀x ∈ Z 2 , g ∈ G) or graphically , (2.10) and the projectors obey that is .
(2.12) Considering again the case of a compact Lie group, the first condition in (2.9) implies j s ⊗ j t ∼ = j l ⊗ j d , i.e. the congruence between combined physical representations and combined ingoing virtual ones.The remaining two conditions imply j s ∼ = j r and j t ∼ = j u , meaning that the individual physical representations are identical to the virtual ones on the same leg.

Transfer operators
The key element for our entanglement calculations in the next section is the so-called transfer operator.The local on-site transfer operator is defined by taking two copies of the physical-virtual tensor and contracting the physical indices between them.This general prescription is valid for any tensor network type.We employ it here to the gauged PEPS framework, for which the local on-site transfer operator can be written as (2.13) This tensor T (1) can be interpreted as a map from the ingoing Hilbert spaces (represented by ket vectors on the l and d legs) to the outgoing ones (given as dual bras on the r and u legs).The matrix elements of T (1) are given as Here, two sets of physical indices, {s, t} and {s ′ , t ′ }, of a physical-virtual PEPS tensor A st ruld and its adjoint Ās ′ t ′ r ′ u ′ l ′ d ′ are contracted, leaving a tensor with in total 8 indices, {ll ′ , rr ′ , dd ′ , uu ′ }.
Let us study properties of the transfer operator when imposing the gauge symmetry.The symmetry constraints (2.9) of the PEPS tensors imply . (2.15) That means that vectors on the ingoing legs combine together to a joint singlet under (θ l g ⊗ θl ′ † g ) ⊗ (θ d g ⊗ θd ′ † g ), while vectors on the outgoing ones are separately on-leg singlets under the respective transformations (θ g ⊗ θ † g ).Symbolically, we can illustrate this map as follows . (2.16) When calculating entanglement quantities in the next section, we tile the whole 2D lattice with transfer operator combinations.From (2.16) it becomes apparent, that in such a case the input legs (l, d) are contracted with on-leg singlets from the right and up index of the neighboring site to the left and below.Defining the projector as Π 0 ≡ j |0(j)⟩ ⟨0 (2.17) As shown in [95], (τ 0 ) j l ,jr;j d ,ju is a positive symmetric matrix, for which a decomposition with orthogonal V and diagonal Λ (containing the eigenvalues λ µ ) is existing.The operator τ0 hence can be written as τ0 = where the Mµ are given as Mµ = j 1 ,j 2 V j 1 j 2 µ |0(j 1 )⟩ ⟨0(j 2 )|.In (2.19), these act independently on the horizontal (l, r) and vertical indices (d, u).
3 Area laws for Rényi entanglement entropies

Idea in 1D
In this section, we introduce the underlying idea of a transfer operator approach, which allows us to calculate normalized Rényi entanglement entropies of order n > 1, defined as Here, the reduced density matrix is defined as ρ A = Tr B [|Ψ⟩ ⟨Ψ|] for a pure quantum state |Ψ⟩.Contrary to the usually defined Rényi entanglement entropy, we assume an additional normalization in (3.1), allowing us to simplify the construction for the subsequent PEPS framework.Note that in the limit n → 1 + , S n approaches the entanglement entropy, defined as the von Neumann entropy of the reduced density matrix.
The only nontrivial element in the definition (3.1) is the trace in the numerator.We demonstrate this calculation for the case n = 2.For a 1D QMB state, which is expressed in MPS form (and without assuming any gauge invariance), the so-called purity takes the form In eq.(3.5), we make use of diagrammatic MPS notation. 5Each row corresponds to the alternating ket-and bra-tensors of the quantum state.We assume an infinitely large, periodic chain and exemplarily defined the subsystem A as the 3 tensors in the middle, such that the complement B is given as the outside region (cf.green marks).Different types of traces in both regions are indicated by connected physical indices of MPS tensors.The crucial observation in this expression is the following.The quantity of interest is uniformly composed out transfer operators.Inside the subsystem A in the tensor network diagram in (3.5), we have marked a single-site transfer operator E (1) by a small red box.In fact, this tensor appears twice, such that one can write an entire column as E (2) ≡ E (1) ⊗ E (1)  (indicated by a large red box).In the complement region B, we have marked a repeating tensor Ẽ(1) , which is marginally different from E (1) by the ordering of indices.
A nontrivial tensor contraction appears only at the two subsystem boundaries, marked by green dashed lines in (3.5).Let us consider exemplarily the right one between subsystem A and the complement B (cf. the blue box).This part of the tensor network diagram can be simplified as follows: .
For the first equality in (3.6) we have rewritten the expression in terms of E (1) and Ẽ(1) .In the next equality we made use of the fact that Ẽ(1) is given as Ẽ(1) = P E (1) P , where P is the operator permuting the two tensor indices on each side of the transfer operator.We finally defined the boundary operator X, which realizes the shifting and permutation of tensor indices in-between two copies of E (2) .If we label the tensor indices from top to bottom as 1, . . ., 4, X can be effectively seen as the following permutation operator One easily finds that the left boundary of the subsystem takes the same form.Since all remaining permutation operators to the left and right side cancel each other, the purity can be re-expressed entirely in terms of E (2) and two insertions of boundary operators X, This expression, and its generalization to higher orders n, allows us to calculate the desired Rényi entanglement entropies in (3.1).

Contraction of transfer operator plaquettes
The previously introduced transfer operator method can be directly generalized to two dimensions using PEPS.We explain this method again for order n = 2 here.As it became apparent from the MPS construction, the essential tensors for calculating the purity are a double copy of the transfer operator and the boundary operator.The structure of the latter is unaltered.For the former, we make use of the operator τ0 , which takes the gauge symmetries into account (cf.section 2.2), and define a single-site double copy of the transfer operator as T (2)  . (3.9) Here, we have introduced a new plaquette notation for T (2) , represented by a gray square.
When calculating the purity and subsequently other entanglement measures for a LGT on a 2D lattice, we have to define T (2) for each lattice site.The resulting tensor network setup is illustrated in Fig. 3.We assume a periodic system of N 1 ×N 2 lattice sites x ≡ (x, y) in the horizontal and vertical direction, respectively.The lattice is tiled with transfer operator plaquettes T (2) .In general, we assume that a subsystem A can be defined as any connected region within this 2D lattice, which includes only full plaquettes.In Fig. 3, we have placed a rectangular subregion A in the lower left corner of the lattice (marked as a green frame).The subsystem has size R 1 × R 2 , where we measure the dimensions R 1,2 in units of the lattice spacing a. 6 Along the boundary of A, the boundary tensor X is inserted at each lattice site.
To calculate p 2 one needs to contract all tensors in both dimensions.This procedure can be done either column-or row-wise.For the following evaluation we choose to calculate transfer rows by contracting tensors in the horizontal direction first.The transfer rows are afterwards contracted in the vertical direction from bottom to top.This yields the following expression for the purity (3.10) In eq.(3.10), we have defined three types of transfer rows, which are visualized in Fig. 3.
In the bottom row, the boundary operator X is inserted at each lattice site 0 < x < R 1 along the boundary of A. We define this operator as where outside the subsystem (x > R 1 ) only identity operators are inserted (indicated by a green dashed line in Fig. 3).The rows above, which are inside the subsystem A, include T (2) at each site and one insertion of X at the left and right boundary.We denote the corresponding transfer operator as E || , taking the form depends on R 1 through the second insertion of X at position x = R 1 .After one additional insertion of X (R 1 ) at the top of the subregion, the tiling of the 2D lattice is completed by the remaining copies of the transfer row E (2) consisting solely out of singlesite transfer operators, Formula (3.10) is derived by taking R 2 /a copies of E || and N 2 − R 2 /a copies of E (2) to tile the whole lattice.
The normalization Tr[ρ A ] 2 in (3.1) can be easily included by observing that Tr 2 .Defining a single-copy transfer row identically to the Wilson loop studies in [95] as we have (3.15)

Hilbert space decomposition and Gauss laws
As alluded in the introduction, our transfer operator approach implicitly assumes an extended Hilbert space construction [77,78,[81][82][83] for the calculation of the normalized Rényi entanglement entropies.In this framework, the tensor product of the gauge field Hilbert spaces H l on all links l forms an embedding for gauge-invariant states.A partial trace over the link Hilbert spaces belonging to the complement region B is well-defined, such that the reduced density matrix can be calculated and the system admits the desired decomposition Apart from the Hilbert space decomposition, it is important to analyze the gauge invariance for this approach, which is encapsulated in the Gauss laws.In the context of entanglement measures, the relevant quantity is the reduced density matrix ρ A .We want to exemplarily discuss gauge invariance of ρ A for the Abelian U (1) group in this section.Naturally, modifications of the Gauss laws can be expected only at the subsystem boundary.Moreover, lattice points in B, which are not neighboring with A, are not relevant for the following discussion.
The precise setting of the link decomposition follows from the plaquette tiling of transfer operators.As defined in the previous section and illustrated in Fig. 3, the subsystem A consists of full plaquettes.A single transfer operator τ 0 and its double copy T (2) in (3.9) contain the gauge fields on the right and up leg.As a consequence, all four gauge fields around any lattice site x ∈ A ⊂ Z 2 , which is not at the left or bottom boundary of the subsystem, are completely inside A. For those lattice points, the reduced density matrix is invariant under where the Gauss law operator (in absence of charges) takes the standard form Here, we have labelled the electric field variables E i as indicated in Fig. 4(a).
The first exception to this rule appears for lattice sites inside A at the left and bottom boundary of the subsystem.This is due to the fact that at these locations, the gauge fields around a lattice site belong to different regions.(Compare Fig. 3, where we have marked the gauge fields along the boundary of A by blue dots.)At these parts of the boundary, three different types of exceptions arise, which are illustrated in Fig. 4(a-c).In case (a) the subsystem boundary is at the left side of the plaquette, hence the third gauge field belongs to region B, while all others are inside A. Similarly, in case (b) for the lower left corner, E 3 and E 4 are outside of A, while for case (c) only E 4 is inside B. The invariance properties of ρ A at these points can be derived by using the gauge invariance of the full density operator ρ = |Ψ⟩ ⟨Ψ|, given as e iϕG(x) ρ e −iϕG(x) = ρ w.r.t. the full Gauss law operator (3.18) (∀x ∈ Z 2 ).In case (a) we therefore have In the second line we have used the fact that only E 3 is inside B, such that the other electric field values can be pulled out of the trace.In the following line we made use of the cyclicity of the trace, canceling the phase factors.The final expression in eq.(3.19) takes the usual form of a gauge invariant operator, but w.r.t. the newly defined effective Gauss law operator Using the same methodology, we immediately deduce that for case (b) the gauge transformation takes the same form with an effective Gauss law operator given as G(x) = E 1 + E 2 , while in case (c) we have The second type of Gauss law modifications appears for lattice sites along the top and right subsystem boundary.Here, only the neighboring points in the complement region x ∈ B are affected.Panel (d) and (e) in Fig. 4 show the two possibilities.The gauge transformation of ρ A is again given by the structural form (3.19).In case (d) (representing the right boundary), G(x) follows as G(x) = −E 3 , and in case (e) (representing the top boundary) one has G(x) = −E 4 .
In summary, we have discussed the principles of gauge invariance in the proposed transfer operator approach, which is based on an extended Hilbert space formalism.We have shown how the gauge invariance of the full density matrix, given by the regular Gauss laws (for all lattice sites), implies modified Gauss laws for the reduced density matrix at lattice sites along the subsystem boundary.The latter can be represented using an effective Gauss law operator.For notational convenience we have chosen the U (1) gauge group in the previous discussion.The framework, however, can be directly generalized to arbitrary other gauge groups.

Results in the thermodynamic limit and continuum
In this section we derive long-range properties of the normalized Rényi entanglement entropies (3.1) using our transfer operator setup.The results are obtained by considering the thermodynamic limit for a large subsystem within an infinite lattice.Analogous to the Wilson loop analyses in [95], we make use of an eigenvalue decomposition of the individual transfer matrices, defined as i ⟩ ⟨w i | (3.21) i ⟩ ⟨w In these expressions we have considered the diagonalization of all transfer matrices E (1) , E (2) , E || , and denoted their eigenvalues respectively as ρ i , ρ i , ρ i .The corresponding right eigenvectors are given by the kets |v i ⟩, while the associated bras (labelled as w) denote the corresponding left eigenvectors.Since E (2) || (R) spatially depends on the insertion of the boundary operator, its eigenvalues and -vectors also depend on R. By construction, the eigenvalues and -vectors of E (2) are related to those of E (1) and |v i ⟩ .
We assume all eigenvalues to be ordered decreasingly, and the existence of a spectral gap with || ), i.e. |ρ Using these diagonal forms of the transfer matrices, we can rewrite the purity (3.10) as In this expression, we have factored out the (possibly degenerate) dominant eigenvalues ρ 1 of the two relevant types of transfer rows.For large loops, as assumed in the thermodynamic limit (3.20), the second sums containing the subdominant eigenvalues in (3.25) vanish.Moreover, eq.(3.25) includes a factor dim 4(R 1 +R 2 )/a (J), which originates from the number of contracted indices along the subsystem boundary, each taking dim(J) values.As argued in [95], the singular values depend only on the irrep J, but not on these indices.Hence all dim(J) are equal and contribute only the overall factor to the result.
Let us now include the normalization factor by defining the normalized purity

Using (3.24), we observe that (Tr
, which cancels the dependence on N 2 in (3.25).Employing all simplifications and taking the overall trace, we are left with To further simplify the resulting equation for p2 , we make the important observation that we alternatively, and fully equivalently, could have tiled our lattice column-wise, instead of the row-wise contraction scheme.Under this assumption, p2 is calculated as (compare with eq.(3.10) Here, X (R 2 ), E || (R 2 ), E (2) and E (1) are calculated as defined in eqs.(3.11), (3.12), (3.13) and (3.14), respectively, but under the replacements N 1 → N 2 , R 1 → R 2 and a column-wise contraction (i.e.trace).Following the same steps that lead to (3.27), we get The result (given as the overall trace for the purity) is invariant under the order of the contraction scheme, hence formulas (3.27) and (3.29) have to be identical.In Appendix A we prove that this condition can be only fulfilled if the dominant eigenvalue ρ 1 (R) has, in the most general case, the spatial dependence ρ ′( 2) where Γ (2) and κ (2) are real coefficients.Upon identification of the two formulas, the normalized purity is shown to follow as Notably, this structural form of p2 is very similar to that of the Wilson loop expectation value found in [95].In particular, the dominant eigenvalue ρ of a single transfer row with two gauge field insertions was found to obey the same structural decay behavior as in (3.30), i.e. ρ ′(1) 1 = Γ (1) e −κ (1) R .Moreover, the distinction between the two possible cases thereof -spatial independence (κ (1) = 0) and exponential decay (κ (1) > 0) -leads to the Wilson loop perimeter law and Wilson loop area law, respectively.In the present context, we are instead interested in the behavior of the normalized Rényi entanglement entropy and want to explore the same corresponding regimes for κ (2) .From the definition (3.1), we have The first term in (3.32) is simply a constant.The last term depends on the perimeter of the subsystem, which is multiplied with the logarithm of the ratio of dominant eigenvalues.As such, it is a manifestation of the entanglement area law.The second term, on the other hand, depends on the area of A and hence would represent an entanglement volume law.
From our physical expectations and the PEPS construction, we can exclude the existence of such a term: For a PEPS with bond dimension χ and a subsystem with perimeter P (number of boundary lattice sites), one can easily show (using eqs.(3.2) and (3.3)) that S n ≤ S 1 ≤ P ln χ.The PEPS ansatz naturally implements the entanglement area law.We therefore can conclude that the dominant eigenvalue of the transfer row E (2) || may not have an exponential decay, i.e. ρ ′( 2) does not depend on the dimensions of the subregion, ρ 1 (R).Using the notation ρ ′(2) 1 ≡ Γ (2) = const, we write the final result for S2 as .
(3.34) Formula (3.34) is the first main result of this paper.It allows to calculate the second normalized Rényi entanglement entropy in the thermodynamic limit for an arbitrary 2D LGT in terms of properties of transfer matrices.In contrast to the complementary results for the Wilson loop expectation value in [95], which are also based on a transfer operator construction, we here cannot distinguish between two different geometric regimes.While the Wilson loop perimeter vs. area law detects the presence of a deconfined or confined phase, the Rényi entropy always exhibits an entanglement area law, and has no other apparent phase differences.On these general grounds, the Rényi entanglement entropy thus cannot serve as an equivalent probe of (de)confinement in pure LGTs.Nevertheless, as we will demonstrate in the next section, (de)confinement properties do leave an imprint on the behavior of entanglement.This, however, requires the study of parameter regimes for specific LGTs.
Some further properties of S2 can be concluded from (3.34).First, since the entropy is manifestly positive, S2 ≥ 0, we infer ρ 1 , which enforces a positive contribution through the logarithmic factor.The second important observation concerns the opposing infrared (IR) and ultraviolet (UV) limits.From our transfer operator derivation we see that an infinitely large subsystem (R 1,2 → ∞) causes an IR divergence in (3.34), while in the continuum limit a → 0 an UV divergence emerges.Both are expected features of entanglement quantities in QFTs [113].

Generalization to higher orders
The transfer operator construction can be directly extended to calculate normalized Rényi entanglement entropies (3.1) of (integer) order n > 2. For that purpose we consider the quantity which is the higher-order generalization of the purity (3.4).For a MPS, p n is structurally constructed in the same way as exemplified in (3.5) by adding the corresponding layers of the quantum state |Ψ⟩ and its adjoint ⟨Ψ|.As for p 2 , the resulting tensor network for p n is uniformly composed out of transfer operators and two insertions of boundary operators X (n) .Repeating the same scheme as in (3.6) for the right boundary, one finds that X (n) acts as the following index map As in section 3.1, we have labelled the indices from top to bottom as 1, 2, . . ., 2n, and assumed periodicity for the last one, i.e. {2n} → {2}.At the left boundary, the mapping is inverted, which is described by the transposed operator X (n)T .By contracting the whole tensor network, p n is calculated as the trace.
We can directly generalize this construction to gauged PEPS by repeating the idea of lattice tiling.Defining the subsystem as before, and including the normalization factor, we then have Here, the basic element is the plaquette operator which contains n tensor product copies of the transfer operator.Based on this construction, we have the following generalized definitions: As in the previous section, we can make use of the eigenvalue decomposition and demand the invariance of the contraction scheme under the order.This directly leads to the formula Expression (3.39) is valid in the thermodynamic limit for any LGT PEPS.It generalizes (3.34) and is the second main result of this paper.As in the previous discussion, the dominant eigenvalue ρ || may not depend on the subsystem dimensions to forbid an entanglement volume law.For the dominant eigenvalue ρ holds.It becomes discernible that the previous conclusions for the entanglement behavior hold equally for arbitrary Sn .That is, eq.(3.39) manifests the entanglement area law for all normalized Rényi entanglement entropies.The results persist in the continuum limit, in which an UV divergence appears when a → 0 is taken.
These general results do not yet allow any conclusions regarding entanglement properties in the deconfined versus confined phase.In particular, we would expect larger amount of entanglement in the deconfined phase.The Wilson loop results in [95] shared that expectation through the observation that the eigenvectors of the transfer matrix E (1) are farther away from being product vectors in the deconfined phase.The result (3.39) would support that hypothesis, if the ratio of dominant eigenvalues, ρ 1 , tends to decrease in the deconfined phase.However, from the general formalism outlined so far, we cannot conclude this behavior.It is therefore necessary to explicitly apply the transfer operator approach to a specific LGT and study its properties in the different phases.
As an alternative to the Rényi entanglement entropies, one could also compute the entropic c-functions [114][115][116], which extract subleading corrections to the area law term of the entropies themselves.This calculation, however, is only meaningful if the entangling surface is independent from the size of the subregion.Such a scenario exists for example when one considers an infinitely long slab as A (with width R 1 and length R 2 → ∞), which differs from our assumption in (3.20).The entropic c-function is then proportional to the derivative of S n w.r.t.R 1 .This quantity is UV-finite and counts the degrees of freedom at the respective length (or energy) scale [116].Such a setup has been implicitly employed in previous holographic explorations, and could also be used to study the limit of large-N c gauge theories in an extension of our present analysis.Interestingly, in the recent work [117], the c-function was calculated for the classical 3D Ising model, which is dual to the Z 2 LGT that we study in the next section.4 Entanglement and (de)confinement for the Z 2 LGT In this section we apply the transfer operator approach to the Z 2 LGT PEPS as an explicit example.By comparing with the Wilson loop results of [95], we connect the entanglement characteristics to (de)confinement properties using numerical calculations.

Transfer operator setup
We consider a Z 2 LGT, for which the physical gauge fields are placed on the links of a 2D lattice.The two-dimensional group Hilbert space on each link is spanned by the spin representations j = +, −.Under the action of the x-Pauli matrix as the Hermitian group element operator, spins are inverted, i.e.Θ(x) |ψ⟩ = |ψ⟩ ∀x ∈ Z 2 .Here, the z-Pauli matrix is the group operator acting as σ z |±⟩ = ± |±⟩.The second group operation is represented by the trivial identity operator.
In contrast to the most general form (2.1) of the gauge transformation, we do not have any difference between left and right group operations here.We consider the most general gauged PEPS ansatz, which is translationally and rotationally invariant.Moreover, we assume a minimal construction with the smallest possible bond dimension, i.e. physical and virtual spaces are identically spanned by the spin states |±⟩.As shown in [95], there are only 4 parameters necessary to fully encompass all nonvanishing gauged PEPS tensors A st ruld (s, t, r, u, l, d = ±) as follows: The single-site transfer operator T (1) is then constructed in the same way as in [95] using (2.13) and (2.14).In fact, we can directly infer the operator τ0 , which takes the simplifying gauge symmetry constraints into account when calculating traces in the 2D lattice, by identifying the vector space elements of the transfer operator legs.
where in this convention the rows correspond to the horizontal PEPS directions (l, r), while the columns correspond to the vertical ones (d, u). 8 By diagonalizing (4.3), we can rewrite τ 0 in the form (2.18).The matrix V takes the form where the entries v 11 , v 12 , v 21 , v 22 are determined numerically for chosen values of the PEPS parameters. 9With the eigenvalues and the operators the transfer operator then takes the desired form (2.19). 10s outlined in the previous section, when calculating purities and Rényi entanglement entropies, we need to construct transfer rows.For the following studies, we restrict ourselves to the second order (n = 2), and calculate the normalized Rényi entanglement entropy S2 .Using the above setup, we can construct the operator E (1) directly as where the indices µ x = 1, 2, 3, 4 are taken for each lattice site x = 1, . . ., N 1 along a row.The operator E (2) = E (1) ⊗E (1) follows directly as a tensor product copy of E (1) .However, for numerically evaluating this expression in matrix form, it is more convenient to choose the tensor product structure site-wise.We therefore have where now the set of both {µ} and {ν} indices run over all values independently.For the transfer row E || , the single-site boundary operator X, given as is additionally inserted at positions x = 1 and x = R 1 as follows (4.9) from the subsystem width R 1 .Right: Dependence of the normalized Rényi entanglement entropy S2 on the subsystem length R 2 for several values R 1 (colored dots).The colored curves represent the linear prediction in the thermodynamic limit for each case, confirming the entanglement area law of the numerical data.See text for further explanations.
Defining the boundary row X (R 1 ) as we can numerically calculate the normalized purity as follows The normalized Rényi entanglement entropy is finally given by S2 = − ln p2 .(4.12)

Numerical results
In full generality, the outlined Z 2 LGT PEPS setup has a 4-dimensional complex parameter space {α, β, γ, δ} ∈ C 4 .The analytical and numerical Wilson loop studies in [95] revealed that it is in fact a nontrivial task to identify the confined phase in this theory through a Wilson loop area law.In particular, it was shown that only in the perturbative regime |β| ≪ |δ| ≲ |α| and γ = 0, one finds a Wilson loop area law, i.e. a confining phase, as long as |δ| ̸ = |α|.In the case |δ| = |α|, there exists a deconfined phase.In the first part of our numerical studies, we want to confirm the existence of an entanglement area law as found in our general analytical treatment.In all subsequent analyses, we restrict ourselves to the real and positive parameter space R 4 + .We consider the confining phase in the perturbative regime and define a thin but long 2D lattice of dimensions N 1 = 4, N 2 = 100, and set α = 1, β = 0.1, γ = 0, δ = 0.95.We implicitly assume unit lattice spacing a = 1.
The left panel in Fig. 5 shows the dominant eigenvalue ρ || (R 1 ) (orange dots) for three different subsystem widths R 1 = 1, 2, 3.The numerical values are identical for all cases and hence confirm our theoretical prediction that ρ As a comparison, we also plot the dominant eigenvalue ρ (2) 1 of E (2) (blue dots), which allows for a positive contribution to S2 in the analytical formula (3.34) 1 .The right panel in Fig. 5 shows the resulting values of S2 for a large range of subsystem lengths R 2 and several parameter values R 1 (colored dots).On top of the data points, linear functions are shown for each parameter, which follow directly from the analytical result (3.34) when the previously determined numerical values of ρ are inserted.The numerical data points agree perfectly with the linear prediction.These results confirm the presence of an entanglement area law.Moreover, it demonstrates that the chosen numerical setup already reproduces results in the thermodynamic limit, i.e. for infinitely large systems.These conclusions are independently also valid in other parameter ranges of the Z 2 LGT PEPS when a nonperturbative or deconfined regime is considered.
Our analytical studies in the previous section have shown, in contrast to the Wilson loop expectation value, that there is no geometric distinction between the confined and deconfined phase for the Rényi entanglement entropies in 2D LGT PEPS.To study the interplay of (de)confinement and entanglement in this situation, it is instead necessary to consider the dependence on the underlying parameter space of the theory.For that purpose we study in Fig. 6 the dependence of S2 on δ for α = 1 and β = 0.1, i.e. in the perturbative regime as before.At γ = 0 (blue curve), S2 is monotonously increasing on a very small scale as long as δ < 1, while it is rapidly increasing for δ > 1.At δ = 1, there is a nonanalytic kink in the data.This point α = δ agrees precisely with the singular deconfined point in this regime.In other words, deconfinement leaves a clear imprint on S2 .However, it does not match the expectation that the deconfined phase should be associated with a larger value of S2 since the amount of entanglement is expected to increase.
One important result of the Wilson loop study in [95] was the following conclusion: In the confining phase at γ = 0 of the perturbative regime, the eigenvectors of the transfer row E (1) are product vectors.If γ is switched on, the eigenvectors are taken further away from product form; the Wilson loop area law is broken and the phase becomes deconfined.
In Fig. 6 we study the corresponding behavior of the Rényi entanglement entropies for such a scenario.The orange and green curve show the behavior of S2 for γ = 0.1 and γ = 1, respectively.While the γ = 0.1 curve is nearly identical to the previously discussed one at γ = 0, S2 deviates for the largest plotted value of γ.Specifically, the nonanalyticity at δ = 1 disappears in the latter case and the values are lower in the regime δ > 1.
In Fig. 7 we extend these analyses by showing S2 in the whole γ − δ parameter plane.The left plot is in the perturbative regime at β = 0.1, and the right one in the nonperturbative regime at β = 1.For the perturbative case it becomes discernible that there is a nontrivial structure with local minima and maxima as γ is increasing (i.e. when the theory becomes deconfined).The previously observed rapid increase of S2 towards large values of δ persists only as long as γ ≲ 1.In the nonperturbative regime, i.e. for a completely deconfined phase, S2 is instead monotonously decreasing towards larger values of γ.The conclusions of these numerical analyses are the following.First, our transfer operator approach allows us to reliably calculate Rényi entanglement entropies for arbitrary LGTs in two dimensions, which we demonstrated for the example of Z 2 .In fact, we could achieve results for infinitely large systems in the thermodynamic limit already for moderate system sizes.Second, when studying the interplay of (de)confinement and entanglement, it is necessary to derive and analyze the PEPS parameter space of the considered LGT.We then found that (de)confinement does leave an imprint on the behavior of the Rényi entanglement entropy for specific parameter dependencies.However, it becomes also clear that these entropies cannot serve as an equally reliably order parameter of (de)confinement as the Wilson loop expectation value.In particular, it is not discernible that S2 (or any higher order entropy) explicitly counts the number of degrees of freedom in our specific setup, which is in contrast to higher-dimensional results for the entanglement entropy S 1 (cf. the discussions in the introduction).

Discussion and outlook
Motivated by recent explorations in quantum information science, particle physics and holography, we have studied the interplay of entanglement and confinement for 2D pure LGTs in this paper.We have developed a transfer operator approach based on PEPS, which allowed us to calculate normalized Rényi entanglement entropies (3.1) within an extended Hilbert space formalism.The basis of our approach is the observation that the relevant tensor network diagram is uniformly composed out of transfer operators and additional operator insertions along the subsystem boundary.The whole lattice therefore can be tiled by plaquettes of transfer and boundary operators.Using the gauge invariance of the full pure state density matrix, we then have analyzed modifications of the Gauss laws for the reduced density matrix along the subsystem boundary for this approach.
Local properties of transfer operators on-site and along rows of the lattice turn out to dictate the long-range properties of the entanglement entropies.This is encapsulated in our final result (3.39), which is valid in the thermodynamic limit.It exemplifies the entanglement area law and exhibits an UV divergence in the continuum limit.The transfer operator approach hence provides a powerful tool to study entanglement properties in arbitrary pure LGTs.However, in contrast to the Wilson loop expectation value as another nonlocal observable, the normalized Rényi entanglement entropies do not have different geometric regimes or other features indicating the presence of a confined or deconfined phase.
We have explicitly applied our methodology to the Z 2 LGT PEPS with minimal bond dimension.Using numerical evaluations, we reproduced results for S2 in the thermodynamic limit already for moderate lattice dimensions.By comparing with Wilson loop results of [95], we found that (de)confinement does leave imprints on the behavior of the Rényi entanglement entropy in specific parameter regimes.However, the entropies are not an unambiguous probe of confinement or deconfinement in the whole parameter space of the theory.
For future studies, it would be highly interesting to include dynamical matter (i.e.fermionic degrees of freedom) to the setup.In this for the standard model relevant physical situation, the structure of the tensors and transfer matrices would be modified.Recent results based on a gauged Gaussian fermionic PEPS ansatz [66], which employs sign-problem free Monte-Carlo techniques for tensor contractions [67], could provide a promising starting point for that purpose.Moreover, the first explorations of this ansatz in 3+1D [118] could be used to generalize our approach to higher dimensions.

Figure 1 .
Figure 1.For a pure (2+1)D LGT in the Hamiltonian formalism, the gauge fields are placed on the links of a spatial lattice (indicated by the blue lines).When acting with the gauge transformation (2.1) on all four links around any lattice site x, the represented quantum state remains invariant, cf.eq.(2.2).

Figure 2 .
Figure 2. Ingredients of a gauged PEPS tensor network.(a) On each lattice site a tensor A st ruld with 4 virtual indices r, u, l, d and two physical gauge indices s, t is placed.(b) Virtual states on neighboring sites are projected onto maximally entangled states via contraction with tensors B 1,2 on the links.(c) A PEPS with only physical indices (after projection) is obtained on the lattice.

Figure 3 .
Figure 3. Setup for the calculation of the purity and Rényi entanglement entropy for a LGT on a 2D lattice.The lattice is tiled with transfer operator plaquettes (gray squares) defined in (3.9).A rectangular subregion A is chosen in the lower left corner (green frame).The purity (3.10) is calculated by contracting different types of transfer rows (orange boxes) and boundary operators (green dashed lines).Blue dots indicate the position of gauge fields at the boundary of A. See text for further explanations.

Figure 4 .
Figure 4. Modifications of the Gauss laws along the subsystem boundary: The transfer operator approach is based on an extended Hilbert space formalism, which allows to realize a Hilbert space decomposition.Panels (a-e) illustrate all different possibilities of lattice sites along the boundary of the subsystem A. In these cases, the reduced density matrix ρ A is gauge transformed according to (3.19) for an position-dependent effective Gauss law operator G(x).Blue dots indicate gauge fields on links around a lattice site x. Green lines mark the subsystem boundary to the complement region B. A gray frame indicates a single transfer operator plaquette (cf.also Fig. 3).See text for further discussions.

Figure 6 .
Figure 6.Dependence of S2 on the parameter δ in the perturbative regime.The data demonstrate that (de)confinement leaves an imprint on entanglement properties: At γ = 0 (blue curve), the theory is confined except at the singular point α = δ = 1 as signaled by a nonanalytic kink.The kink disappears towards larger values of γ (orange and green curve) when the theory becomes deconfined in the whole parameter space.See text for further explanations.

α 1 Figure 7 .
Figure 7. Behavior of S2 in the γ − δ parameter plane.The left plot is in the perturbative regime, the right one in the nonperturbative regime.See text for detailed discussions.