Wilson Loops and Area Laws in Lattice Gauge Theory Tensor Networks

Tensor network states have been a very prominent tool for the study of quantum many-body physics, thanks to their physically relevant entanglement properties and their ability to encode symmetries. In the last few years, the formalism has been extended and applied to theories with local symmetries to - lattice gauge theories. In the contraction of tensor network states as well as correlation functions of physical observables with respect to them, one uses the so-called transfer operator, whose local properties dictate the long-range behaviour of the state. In this work we study transfer operators of tensor network states (in particular, PEPS - projected entangled pair states) in the context of lattice gauge theories, and consider the implications of the local symmetry on their structure and properties. We focus on the Wilson loop - a nonlocal, gauge-invariant observable which is central to pure gauge theories, whose long range decay behaviour probes the confinement or deconfinement of static charges. Using the symmetry, we show how to handle its contraction, and formulate conditions relating local properties to its decay fashion.


I. INTRODUCTION
In recent years, tensor network states [1] have been a very prominent tool, rooted in quantum information science, for the study of quantum many body systems and especially strongly correlated physics. In particular, Matrix Product States (MPS) [2,3] enable to study numerically and analytically physically relevant states, e.g. ground states of local many body Hamiltonians (that is, states exhibiting an entanglement entropy area law). In higher spatial dimensions, MPS generalize to PEPS -Projected Entangled Pair States [3,4]. These are useful for the description of strongly correlated physics in two spatial dimensions and more.
PEPS (and MPS) are constructed out of the contraction of local building blocks (tensors). They satisfy, by construction, the entanglement entropy area law (focusing on the physically relevant part of the Hilbert space) and allows the state to depend on very few local parameters and hence making it feasible for computations (compared with arbitrary states in the exponentially large Hilbert space). Furthermore, it also allows one to encode symmetries already on the single tensor level. By properly parametrizing the local tensors, a global symmetry of the whole PEPS under a symmetry group can be imposed [5,6]. This way, one can generate families of ansatz states in which the symmetry group of the studied model is encoded by construction.
While originally used mostly in the context of condensed matter physics, MPS and PEPS have recently been extended to the study of particle physics too -in particular, to lattice gauge theories, aiming at solving long standing open, non-perturbative questions of the standard model, such as the confinement of quarks [7]. Due to its running coupling [8] Quantum Chromodynamics (QCD) which allows one to use perturbation theory in high energy scales (collider physics) thanks to asymptotic freedom, is strongly interacting in low energies, prevent-ing the use of perturbative methods. Lattice Gauge Theories (LGTs) [7] have been introduced to overcome this difficulty, first as tools for lattice regularization of gauge invariant field theories. They quickly became a very successful numerical approach. Combined with quantum Monte Carlo it has been applied to nonperturbative QCD computations, such as the hadronic spectrum [9]. However, quantum Monte Carlo does not allow for the direct observation of real time dynamics, and faces the fermionic sign problem in several important physical scenarios, not allowing one to probe some of the interesting exotic regions of the QCD phase diagram [10], and this requires the use of other methods, with tensor networks being one such approach. The tensor network framework for lattice gauge theories has been rapidly growing in the last few years.
For 1 + 1d systems, MPS have already been extensively used. This does not only include abstract formalistic descriptions of MPS with a local symmetry [11,12] or benchmarks of models that can be treated in other ways, such as, but not only, the Schwinger model [13][14][15][16][17][18][19][20]. Successful numerical studies of lattice gauge theories in 1 + 1d have been carried out even in scenarios which face the sign problem when approached with conventional methods (such as real time evolution [21][22][23] and finite density [24][25][26][27]). This was done for both Abelian and non-Abelian models -see [28] and references therein for a discussion of that.
The application of tensor networks to higher dimensional lattice gauge theories has been discussed as well in the last few years [29][30][31][32][33][34][35][36][37][38]. From the rather more abstract, or formalistic point of view, gauging mechanisms which lift globally invariant PEPS to locally invariant ones by adding a gauge field and entangling it to the matter properly were introduced and discussed [39,40]. For a parallel approach in the action formalism -tensor field theory -which uses tensor networks (but not tensor network states), see [41] and references therein.

arXiv:2101.05289v3 [quant-ph] 17 Dec 2021
In this work, we will focus on a particular gauging mechanism -the one introduced in [40] and used mostly with fermionic matter, for creating gauged Gaussian fermionic PEPS [31,32]: special PEPS constructions which allow for the description of fermionic matter coupled to dynamical gauge fields. Their construction may be seen as a minimal coupling procedure on the level of states, which is not possible in general but could be done in the context of PEPS [42]. While in general numerical computations are hard and challenging for PEPS in two spatial dimensions and more, it has been shown that, when this particular construction is used, the PEPS may be contracted efficiently (allowing one to extract physical information) when combining with Monte-Carlo methods which do not suffer from the sign problem [33]. Variational Monte-Carlo then allows to find ground states of lattice gauge theory Hamiltonians when such states are used as ansatz states, which has already been demonstrated and benchmarked for a pure Z 3 lattice gauge theory in 2 + 1 dimensions [37].
A question that has to be asked when a PEPS is studied, is how physical information can be extracted from the contracted state -computation of expectation values of observables and correlation functions. Thanks to the special structure of MPS, one may introduce a mathematical object called transfer matrix (or operator) [43] to compute efficiently expectation values of observables and correlation functions. This may be extended to two dimensional PEPS, by first contracting the rows, obtaining effectively a chain of the rows which is an MPS, and considering its transfer matrix [44]. In this work, we will study such transfer operators of lattice gauge theory PEPS in two space dimensions.
Gauge theories are special in the sense that they exhibit a local symmetry, responsible to mediating local interactions between the matter fields. This symmetry gives rise to many local constraints. All the physical states are invariant under gauge transformations -local transformations parametrized by the elements of the socalled gauge group. As a result, only gauge invariant observables and correlation functions -those which are invariant under local transformations -give rise to a nonvanishing expectation values. Thus, LGT PEPS admit a very special structure manifested in the local tensors [30,39,40]; in this work, we focus on the implications of the local symmetry on the transfer operators, and hence aim at using the symmetry to simplify the PEPS contraction, focusing on pure gauge theories (that is, without dynamical matter).
In such scenarios, closed flux loops -usually referred to as the operators which create them, Wilson loops [7] are perhaps the most important observables (and almost the only possible gauge invariant one). The decay rule of large Wilson loops in pure gauge theories serves as probes for confinement of static charges: area law decay implies confinement, while a perimeter law -a deconfined (Coulomb) phase. Confinement implies a gapped, disordered phase, while deconfined phases are massless and ordered [45]. The local ingredients of the Wilson loop are not gauge invariant -only their combination along the nonlocal path preserves the symmetry. This means, that when computing it for a gauge invariant PEPS, the transfer operator formalism must be extended and modified, requiring the inclusion of various types of transfer matrices which construct this nonlocal observable. The different building blocks will also have special properties [32], dictated by the special local symmetry, which will affect the behaviour of the Wilson loop and its decay.
In this work, we will study the properties of transfer operators of gauge invariant PEPS. We will see how the symmetry affects the properties of the local tensors, and that thanks to it, some parts of the tensors may be excluded and ignored when a contraction is done (e.g. when combined with some numerical methods). We will also see how that affects the Wilson loop's decay -that is, how local properties of the tensors dictate the decay of large Wilson loops.
Note that PEPS have been previously used for the computation of Wilson loop expectation value in various cases -Z 2 string nets [46,47], as well as U (1) [31] and SU (2) [32] toy models; here we derive a general framework based on transfer matrix arguments and demonstrate with particular constructions We begin with briefly reviewing important preliminaries from group theory and lattice gauge theory, in section II; move on to formulating gauge invariant PEPS and reviewing their symmetry properties, in section III; in section IV we introduce the transfer operators -after a brief review of their general properties, we formulate the flux-free transfer operators for LGT PEPS, study their properties and use them to calculate the norm; section V focuses on the contraction of Wilson loop expectation values for LGT PEPS, studying the relevant transfer operators and deriving conditions for area and perimeter decay laws; finally, in section VI, we give an explicit illustration, including both analytical and numerical arguments, for a Z 2 lattice gauge theory.
Throughout this work the Einstein summation convention (on doubly repeated indices) is assumed unless stated otherwise; with the only exception of irreducible representation indices, whose summation should not be assumed.
cal theory, we would like to consider transformations parametrized by the elements of G. To do that, for each g ∈ G we introduce a unitary operator θ g , and define it by its action on a basis states of the form |jm . j labels the irreducible representations of G and m is an index labelling all states within this representation -that is, all the states that may be mixed by the transformations θ g that act block-diagonally on the irreps: We can hence express θ g as -and therefore the dimension of the irrep j, dim (j), is also the dimension of H j , the Hilbert subspace spanned by |jm , which we call a multiplet. The Hilbert space may be seen as a direct sum of multiplet subspaces In general, quantum Hilbert spaces may contain more than one multiplet carrying the same irreducible representation. These transformations are sometimes referred to as right transformations, since they mix the multiplet elements |jm , when seen as the components of a dim (j) dimensional vector, via right matrix multiplication, as shown in (1). One can also introduce the left transformationsθ Note that the left transformations are not independent from the right ones: for each g ∈ G one may find h such thatθ g = θ h . We introduce the left transformations separately nevertheless since they will be mathematically convenient later when the PEPS are constructed. When G is a compact Lie group, its elements may be uniquely identified in terms of group parameters or coordinates φ a ; then, for each irrep j, -the parameters φ a (g) depend on the group element, while the generators T j a depend on the representation. The latter form a set of matrices with dimension dim (j), satisfying the group's Lie algebra where f abc are the group's structure constants. One may also introduce the abstract generators, J a , which are block diagonal in the representations, satisfying the algebra too. The states |jm are eigenstates of mutually commuting operators: the j quantum numbers(s) labelling the irreducible representation (and hence the multiplet) are eigenvalues of the Casimir operators which commute with all the generators; within the representation, the states are labelled by the eigenvalues of a maximal set of mutually commuting generators (Cartan subalgebra) -m. Similarly, when the group is finite, j labels the irreducible representation while the m numbers are obtained from the simultaneous diagonalization of a maximal set of commuting transformations.
All the irreps of Abelian groups are one dimensional and thus no m indices are required. In the Z N case, the N different irreps are labelled by the integers j = 0, ..., N −1, which label the group elements g = 0, ..., N −1 too, with D j (g) = exp i 2π N jg . In the U (1) case the group elements are labelled by one parameter as well, φ ∈ [0, 2π), the representations are labelled by integers j ∈ Z, and D j (φ) = exp (ijφ), and T j = j1.
As a non-Abelian example, consider SU (2), whose irreps are labelled by j that are non-negative integers and half-integers. The dimension of each representation is dim (j) = 2j + 1, and the 2j + 1 within the multiplet are labelled by m = −j, ..., j. There are three generators, satisfying the Lie algebra with f abc = abc -the anti-symmetric (Levi-Civita) symbol with a, b, c = 1, 2, 3. The generators in this case are sometimes called the spin or angular momentum components, and then a, b, c = x, y, z. The j = 0 (trivial) representation is one dimensional, with the singlet state |00 . The next representation, j = 1/2, is two dimensional (m = ±1/2 ), with generators proportional to the Pauli matrices, T j=1/2 a the context of quantum simulation and tensor networks, including in this work. Since we consider Hamiltonian lattice gauge theory in 2 + 1 dimensions, our lattice will be two dimensional. As this work focuses on the pure gauge case and matter fields are absent, all the degrees of freedom will reside on the links. We will review the basic ingredients of such models following the conventions of [49] and [50].

Local Hilbert Spaces
Consider a two dimensional lattice, whose sites are labelled by vectors of integers x ∈ Z 2 .ê i denote the unit vectors pointing in directions i = 1, 2, and any link is classified by two numbers, (x, i), standing for the beginning of the link and the direction to which it emanates, respectively. Each link (x, i) hosts a local gauge field Hilbert space H gauge (x, i), which can be spanned by group element states {|g } g∈G labelled by the gauge group elements. These states form a basis of H gauge , with the orthogonality relation where δ (g , g) is the Kronecker delta if G is finite, and a Dirac delta distribution in the compact lie case -denoting the Haar measure of G by dg, Unlike in the multiplet case, here the right and left transformations are independent of one another, as group multiplications: we introduce two sets of unitary operators, Θ g andΘ g , parametrized by the elements of the gauge group G, which implement right and left group multiplications (respectively) on the group element states, The space H gauge can also be spanned by the dual representation basis, whose states are labelled by |jmn -j is an irrep and m, n are identifiers within it. In a sense, using the multiplet states introduced previously, |jmn = |jm ⊗ |jn (13) or, where H j is the dim (j) dimensional subspace spanned by the |jm multiplet states. We read this equation as a decomposition of the link's Hilbert space into a direct sum of products of multiplets of the groups on the left and right of the link, sharing the same irrep. Here one copy of each irreducible representation is used at most; one in the full, Kogut-Susskind case [48], but it is also possible to choose (for example, for reasons of feasibility of computation or experimental implementation) to truncate the sum and not include all the irreps in several ways [49,51] as we will discuss later.
In the non truncated case, using the Peter-Weyl theorem and the group's Fourier transform [49], the transition between the two bases is given by where |G| is the group's volume. In the representation basis, In the compact Lie group case, one can introduce two sets of transformation generators, left and right -L a and R a respectively, such that satisfying the algebra Note that if the group is Abelian, there is no difference between left and right operations and the indices m, n do not exist. Therefore, there R = L ≡ E. Thus, in the U (1) case, for example, we have group states labelled by the single compact parameter |φ and representation states labelled by the single integer |j , related through the Fourier series formula and the representation states |j satisfy For Z N , similarly, we obtain the discrete Fourier series formula While in the SU (2) case, since the group is non Abelian, the situation is more complicated. There are (2j + 1) 2 |jmn states for each j -e.g. one singlet state |000 for j = 0, and four j = 1/2 states, 1 2 , ± 1 2 , ± 1 2 . The group is parametrized by the three Euler angles α, β, γ, and α, β, γ|jmn = √ 2j + 1 8π 2 D j mn (α, β, γ) The Hilbert space in this case is that of a rigid rotator [48,52]. The right and left operators R a and L a correspond to the generators of its rotations in the space and body frames of reference. These two sets of generators commute, and give rise to the same total angular momentum (eigenvalue of the Casimir operator) since it is a rotation scalar quantity which does not depend on the frame of reference [53,54]. Therefore,

Local Gauge Invariance
At each site x, and for each group element g ∈ G, we introduce the gauge transformation which transforms all the four links intersecting at x with respect to the same group element -the outgoing links with the left transformation, and the ingoing ones with the inverse right one. The outgoing links, whose beginning (left) side connects to x, undergo a left rotation, while the ingoing ones, connected through their end (right) side to x, undergo an inverse right rotation.
A gauge invariant state |ψ satisfieŝ (see Fig. 1) and similarly, for a gauge invariant operator O,Θ (one can extend it to the case of static charges [52] which we do not discuss here). In a lattice gauge theory, only gauge invariant states and operators are considered physical.
If G is a compact Lie group, we can formulate the gauge transformationsΘ g (x) in terms of their generators, Gauge invariance is then formulated in terms of the Gauss laws (once again, excluding static charges [52]). We call this eigenvalue equation the Gauss law, since G a (x) can clearly be seen as the divergence of electric fields -L a and R a -on a site. For physical states -the one which satisfy the local constraints (25) -the divergence of electric fields is zero. It is very apparent in the U (1) case, where it takes the explicit form In non Abelian cases, the divergence involves left and right electric fields, which is related to the charge carried by non-Abelian gauge bosons [48] (e.g., the colour charged gluon vs the electric neutral photon).
Since we deal with gauge invariant states, it is expected that the expectation values of non gauge invariant operators will vanish. Thus, when classifying the phases and behaviour of gauge theories one needs to consider only gauge invariant observables and correlation functions.
One option, for compact Lie group, is to compute expectation values of electric field operators and functions thereof (and only of Casimir operators if the group is non-Abelian). Another possible gauge invariant observable is the loop variable, and in particular Wilson Loops [7].
On the local link Hilbert spaces we introduce the group element operators: U j is a matrix of dimension dim (j) × dim (j), whose elements are operators acting on the link's gauge field Hilbert space, H gauge (on each link we can define such operators U j mn ( )). Even though they are Hilbert space operators, all the elements of U j commute -one can see in the definition above that they are all diagonal in the same basis. The matrix elements of U j mix with respect to the transformation properties of the j representation, and, in the compact Lie case, R a , U j mn = U j mn T j a n n L a , U j mn = T j a mm U j m n (32) Let us take some closed path C on the lattice. We define the Wilson loop operator W (C) as the ordered contraction of group element operators along this closed path, that is (33) It is a trace over the product of the group element operators U , seen as matrices, ordered along the closed path C with length L (which is simply the number of links along the path). Depending on the orientation of the path, one may have to use U † instead of U , on half of the links along the path -those pointing leftwards or downwards (see Fig. 2). For simplicity, we will omit the j indices below, but obviously the same irrep must be used along the path, otherwise the matrix product is ill defined. Consider U (1) with j = 1 as an example; there, with half of the phases with a minus sign, according to their orientation. In order to consider the action of the group element operators on representation states, we use the Clebsch-Gordan series and coefficients JM jm|KN [55] and obtain -that is, the action of the group element operator U j on a state with representation J yields states with all representations which are obtained by combining j and J (more precisely, fusing the two irreps together). Acting with a loop operator hence excites the representations along the loop with respect to that rule. One may truncate the Hilbert space in the representation basis: as long as all the irreducible representations used are taken completely and connected by nonzero Clebsch Gordan coefficients when j is added, one may use (35) to define a U j operator acting on that truncated space. The transformation properties (16), (31) and (32) will still hold [49], which may make it convenient for some numerical approaches (or quantum simulation implementations [56]) but, since the group structure will be lost, the group element basis will no longer be defined, making, in particular, (30) and the Fourier transform (15) invalid. In most cases, rectangular Wilson loops are considered. We denote by W (R 1 , R 2 ) a rectangular loop sized R 1 × R 2 (see Fig. 2). Very large Wilson loops of pure gauge theories are a probe for confinement (or deconfinement) of static charges, as introduced by Wilson in [7] (see also [45,57,58]). In a confining phase, for R 1 , R 2 1 (area law), while in a deconfined phase for R 1 , R 2 1 (perimeter law). In [59], Creutz introduced the parameter for the detection of static charge confinement. In the general case of For large R 1 , R 2 , the area factor κ A (called the string tension), should it exist, is the most dominant one. The Creutz parameter χ filters out the contributions of the constant prefactor W 0 and the perimeter coefficient κ P , and thus within a confining phase, χ (R 1 , R 2 ) → κ A > 0 for R 1 , R 2 1, while in a deconfining one it converges to zero.

III. GAUGE INVARIANT PEPS
In this work, we use the lattice gauge theory PEPS formalism of [33,40], with slightly different notations (and restricted to the pure gauge case). First of all, let us review it.

A. Review of the PEPS construction
Each site x ∈ Z 2 of our square, periodic lattice is at the intersection of four legs. The outgoing ones are in the right and up directions, while the left and down directed legs are considered ingoing. We wish, as usual with PEPS, to construct a physical lattice state describing different physical degrees of freedom located on different sites. Each such degree of freedom is described by a local physical Hilbert space: if we had matter, we would fix a physical matter Hilbert space to each lattice site. Here, however, the gauge fields are our only physical degrees of freedom, and they reside on the links. Thus, with each lattice site x we associate two physical Hilbert spaces, located on the outgoing legs. We refer to them as the side (H s ) and top (H t ) physical Hilbert spaces.
These are local gauge field Hilbert spaces (note that the word local here has to do with being defined on a sin-gle link, not with the gauge symmetry being local) -that is, either the full H gauge spaces introduced in Eq. (14), or truncated versions thereof containing only some representations. When truncating, it is important to make sure that all the |jmn for an included j are present, otherwise no gauge invariance can be imposed, as explained above [40,49].
When constructing a PEPS, in order to connect the local physical building blocks to one physical quantum state, one has to introduce auxiliary or virtual degrees of freedom, on top of the physical ones given by the model we study. These are used merely for the purpose of contraction. On each of the four legs we introduce an auxiliary or virtual Hilbert space, H r , H u , H l , H d for the right, up, left and down going legs, respectively. They are spanned by group multiplet states of the form |jm , as defined in Eq. (3). One may include all such multiplets, truncate, or include several copies of the same multiplet, which allows to increase the number of variational parameters; but once again, all the states within a multiplet included must be present, and the representations used in the physical spaces must be included (though possibly with a higher multiplicity). For more details about that, refer to [40] where the general construction of such states is discussed.
On each site, we construct the physical-virtual state where the first ket refers to the physical states and the second to the virtual ones (see Fig. 3). The coordinate x was omitted for simplicity, but the Hilbert spaces are all associated with particular sites and, in general, the tensors A jsmsns;jtmtnt jrmr;jumu;j l m l ;j d m d may depend on the position, although we will focus on translationally invariant PEPS and thus they will be independent of x.
To contract the PEPS, on each link we introduce the maximally entangled states As usual, we construct our PEPS |ψ by projecting the virtual states on the legs onto the maximally entangled states, Note that  of this projection, |ψ , is still a quantum state, including only physical degrees of freedom -the virtual, or auxiliary ones, are all contracted: this is the standard way to contract PEPS, and a frequently used notation. The physical degrees of freedom are now correlated, and in particular, thanks to maximally entangling nearest neighbours, this guarantees the entanglement entropy area law.
One still has some freedom to choose which maximally entangled states to use; the ones that we picked here are invariant under the following group transformations: (with θ g ,θ g defined in (1), (4) respectively) as depicted in Fig. 4. This allows to construct states with a global or local symmetry, as we shall now see.

B. Imposing the local symmetry
We want our PEPS |ψ (42) to be gauge invariant as in (25) with respect to the local gauge transformations defined in (24). If the local physical-virtual states on each site satisfy ( where the physical Hilbert spaces are transformed using Θ,Θ defined in (16), and the virtual ones using Θ,Θ defined in (1) and (4) respectively; see Fig. 5. Using the transformation properties of the maximally entangled states (43) one obtains that |ψ is gauge invariant. In order to get a more intuitive picture of the symmetry conditions (44), let us consider the compact Lie group case again. Omitting the coordinate, since we deal we a single coordinate x, let us denote the right and left generators of the physical degrees of freedom by R s/t a and L s/t a . For the virtual degrees of freedom we can also define such operators, but in their case note that they do not commute, since they do not act on separate degrees of freedom (|jm states, unlike the physical |jmn states). The conditions (44) can be expressed, using these notations, as Gauss laws: The first condition looks like the familiar physical Gauss law. It implies that the two ingoing representations of the virtual indices must combine to the same representation to which the two physical representations combine: j s ⊗ j t ∼ j l ⊗ j d . Therefore, the tensor A jsmsns;jtmtnt jrmr;jumu;j l m l ;j d m d should be proportional to the appropriate Clebsch-Gordan coefficients, The other two conditions, are different identifying the right constituents of the physical degrees of freedom with the virtual states on the same legs. This implies that j r = j s , j u = j t , m r = n s and m u = n t ; A jsmsns;jtmtnt jrmr;jumu;j l m l ;j d m d must be proportional to δ jsjr δ jtjt δ ns,mr δ nt,mu . Combining the first condition with the other two, we can obtain a condition on the four virtual legs: the elements of A jsmsns;jtmtnt jrmr;jumu;j l m l ;j d m d must vanish, unless Examples for constructions satisfying that have been previously given [31,32,40,42]; let us just briefly comment on some special cases. When the group is Abelian, only the irrep indices remain and the Clebsch-Gordan coefficients are simply Kronecker deltas. One can then formulate j r ⊗j u ∼ j l ⊗j d in a very simple way. For U (1), for example, j 1 j 2 |J = δ j1+j2,J , and the Z N is the appropriate modular modification, j 1 j 2 |J = δ j1+j2,JmodN . We thus obtain, in the U (1) case, only tensor elements for which j r + j u − j l − j d = 0 may be nonzero The same applies to non-Abelian groups as well, but since physical states contain the (generally different) m, n quantum numbers it is less simple. For SU (2), e.g., if we choose to include only the j = 0, 1/2 representations, the only non-vanishing tensor elements will be those with an even number of virtual legs (ingoing or outgoing) with j = 1/2, such that a singlet can be formed by combining the contributions of all four legs.
The only freedom left in the definition of A jsmsns;jtmtnt jrmr;jumu;j l m l ;j d m d is to introduce some parameters f j jr,ju,j l ,j d which only depend on the representations, and we obtain [40]: In the following, we will focus on PEPS satisfying the above symmetry properties, with no more than one copy of each irrep in the virtual spaces. This may seem restrictive when attempting to apply the states to real, physically relevant Hamiltonians; here, however, we wish to consider the most minimal constructions which capture the relevant symmetry properties, allowing us to demonstrate our claims and results as accurately as possible. When applied to Hamiltonians as variational ansatz states the states may have to be generalized indeed but in a straight forward way that does not affect the properties we discuss here. For example, as was demonstrated already in the Z 3 case [37], several copies of the virtual representations are required in order to use such PEPS in order to variationally find the ground states of the Z 3 Hamiltonian.
One could also consider a more general PEPS construction, in which such properties are only satisfied after blocking, for effective sites and effective links. The symmetry conditions described above will hold in this case too -for the blocked tensor network, rather than the original, microscopic" one, and thus what we study here could easily be applied to such cases too. A more general scenario would be with local MPO symmetries [4], but this is out of the scope of this work and requires its own, separate discussion.

C. Tensor notation
The projection (42) which generates the PEPS |ψ can simply be seen as a set of contraction rules for the virtual indices of the tensors A jsmsns;jtmtnt jrmr;jumu;j l m l ;j d m d , associating the indices of r at x with those of l at x +ê 1 , as well as u at x with d at x +ê 2 . Hence instead of looking at the local states |A and their projection onto the link stats |B i (x) , we may use, as our basic local building block, where, for the sake of notation simplicity, In (48) the virtual states and their projection are replaced by the matrix products of |l r| along horizontal lines (with the positive direction from the left to the right) and |d u| on the vertical lines (positive direction -upwards). This sets the contraction rules of the tensors A s1s2;t1t2 ruld .
To illustrate, let us reduce to one space dimension and one dimensional PEPS -an MPS [43]. Each local tensor along the one dimensional system includes one physical leg, spanned by states |p , and two virtual ones, on the left and right direction. The state is thus parametrized by the tensors A p lr , and their contraction is simply a matrix multiplication of the virtual indices along the system. For a periodic system with N sites (the modification for open boundaries is straightforward) the state takes the form The PEPS contraction rules in two space dimensions are simply a two dimensional generalization of the trace contraction in the one dimensional case.
The symmetry conditions (44) may also be expressed as properties of the tensor A jsmsns;jtmtnt jrmr;jumu;j l m l ;j d m d . For that, we introduce the (reducible) representation matrices D (g) which are direct sums of the irreducible unitaries D j (g); using them, the symmetry condition (44) may be reformulated as

IV. TRANSFER OPERATORS AND NORMS OF PEPS
Before turning to the study of the transfer operator of our gauge invariant PEPS, let us recall what the transfer operator of a PEPS is. First, we briefly review the one dimensional, MPS case [43]. We strictly focus on the translationally invariant case, since this work is aimed at translational invariant systems; however, the general transfer matrix discussion may be (and has been) generalized to the non-translationally invariant case.
For the computation of an expectation value of some operator O at site x, we will need to define using which we may write Suppose we wish to compute the two-point correlator of O 1 (x 1 ) and O 2 (x 2 ) (assuming for simplicity that We introduce the left and right eigenvectors of E, w i | E = w i | ρ i and E |v i = ρ i |v i , sharing the same eigenvalues ρ i and satisfying the orthonormality relation w i |v j = δ ij , and expand E as Let us sort the eigenvalues in a descending order and assume that the largest one is non-degenerate, that is -the correlations decay exponentially, with a finite correlation length ξ = −1/ log ρ2 ρ1 .

B. Transfer operators of PEPS
The transfer matrix approach can be generalized to two dimensional PEPS, such as the ones we consider here, constructed in (42). We assume the system has periodic boundary conditions -a torus of size N × N (generalizations to other boundary conditions are straightforward). The local transfer operator of a PEPS on a site is a map from two double virtual Hilbert spaces, associated with the ingoing (left and down) legs, to other two double virtual spaces, directed to the outgoing directions (right and up):T = T ll ,rr ,dd ,uu |ll rr | ⊗ |dd uu | (note that we use again a convention in which the input vectors are denoted by bras, in accordance with matrix product ordered from left to right in the positive system directions).
In full analogy with the one dimensional case, the elements of the transfer tensor T ll ,rr ,dd ,uu are given by The norm may be computed by properly contracting products ofT on all the lattice sites; in expectation values of observables, the numerator may be computed by replacingT at the relevant sites by To compute correlations, we will first contract all the tensors along one dimension of the PEPS, converting it effectively to an MPS [31,32,44] whose transfer matrix can be defined as above. For example, the transfer matrix of a row of length N is obtained by contraction along the horizontal direction, Using E and similar transfer matrices which include observables, one may use the entire MPS machinery for computations of norms, expectation values and correlation functions. Naively, one may deduce that correlations in this case decay exponentially as in the MPS case [43]. However, unlike in the one dimensional, MPS case, here the transfer matrix is a composite object with some internal structure, which can lead to different results. It was shown in [60], for example, that two dimensional PEPS can describe critical physics, exhibiting power law contributions.

C. Flux free transfer operators
Let us apply the above to the computation of the norm. For that, consider the flux-free transfer operator, that is, the local building block of the transfer matrix on a single site, with no string (group element operator U j ),T , as defined in (58). We calculate its elements using (59), and thanks to the symmetry conditions (50) we obtain that for every g ∈ G, θ l g ⊗θ †l g ⊗ θ d g ⊗θ †d g T =T θ r g ⊗θ †r g =T θ u g ⊗θ †u g =T (63) (see Fig. 6(a)). This implies, that in (58), the outgoing vectors |rr and |uu are both separately singlets under the action of θ g ⊗θ † g -that is, they are on-leg singlets , denoted by 0 (j r )| and 0 (j u )| and defined as The ingoing legs |ll ⊗ |dd , on the other hand, combine together to a singlet under θ l g ⊗θ †l We can therefore conclude that the general structure ofT iŝ -it is a map with two inputs and two outputs, which takes a joint singlet (on both the ingoing legs) into two separate on-leg singlets, on each outgoing leg alone (see Fig. 6(b)).

D. The row transfer matrix and the norm
Suppose we wish to compute the norm, which involves contracting the tensor product ofT everywhere. EachT obtains its inputs from the neighbour-ingT operators on its left and bottom, whose outputs are on-leg singlets: that is, when the norm is computed, the inputs |j l m l , j l m l 0 (j r )| on the left leg and |j d m d , j d m d 0 (j r )| on the lower one are being contracted with the outputs from neighbouring sites -0 (j r )| and 0 (j u )| respectively. Thus, for the norm contraction it is enough to focus only on a subset of the T elements, where only on-leg singlets are allowed as input. Denoting by Π 0 = j |0 (j) 0 (j)| the projection operator onto on-leg singlets |0 (j) = |jmjm , we definê (introducing a new notation which will be used for tiling diagrams below, in which the legs are implicit). It takes the simple form To see how this simplifies the contraction, let us consider some illustrative examples. First, consider the Z N case, in which (disregarding multiplicities) the virtual Hilbert spaces are spanned by D = N basis states, corresponding to the j = 0, ..., N − 1 irreps. Thus we will have N on-leg singlets, of the form The tensorτ 0 will thus contain N 4 elements; having considered T without taking the symmetry into account, with two N dimensional legs per direction, we would have instead N 8 tensor elements! That is, the number of elements that actually need to be used for contraction is N 4 times smaller. Next, generalize to U (1), and suppose we truncate and allow for the |j| ≤ J for some J > 0. Then we will have once again on-leg singlets of the form (68). There are D = 2J + 1 irreps in the virtual Hilbert space, we have D on-leg singlets and, similarly to the Z N case, we obtain a reduction of D 4 : D 4 elements inτ 0 which we need for the contraction, rather than the D 8 in the most general case.
The simplification is even bigger when we consider non-Abelian groups, because the tensorsτ 0 only see the representations and not the different m values within them. For example, consider SU (2), with the smallest truncation, containing the j = 0, 1/2 representations. This implies that each virtual Hilbert space has dimension 3. Naively speaking, T would be a tensor with 3 8 = 6561 elements. Reducing toτ 0 , with only two on-leg singlets for the two irreps used, the number of relevant elements decreases to 2 4 = 16, that is, approximately 410 times less! If we wish to consider, a little more generally, all the irreps of SU (2) between 0 to some J, the dimension of the virtual Hilbert spaces would be D (J) = J j=0 (2j + 1) = (J + 1) (2J + 1) (note that the sum runs on both integer and half-integer values). Thus T has D 8 (J) = (J + 1) 8 (2J + 1) 8 elements. However the number of on-leg singlets is as the number of irreps, 2J + 1, and henceτ 0 is a tensor with (2J + 1) 4 elements: the reduction factor is (J + 1) 8 (2J + 1) 4 , which scales as J −12 for large cutoffs -a very significant reduction! To examine further the properties ofτ 0 , let us consider (τ 0 ) j l ,jr;j d ,ju as a matrix with the multivalued indices j l , j r and j u , j d . If we assume horizontal-vertical reflection symmetry, we find that it is a symmetric matrix, (τ 0 ) j l ,jr;j d ,ju = (τ 0 ) j d ,ju;j l ,jr Furthermore, it is a real matrix, since using (59) Therefore, there exists an orthogonal matrix V , such that where Λ is a diagonal matrix with eigenvalues λ µ . This allows us to bringτ 0 to the convenient form -one copy of which acts on the horizontal direction and the other on the virtual one. The real matrices M µ form an orthonormal set with respect to the trace inner product. Since V is orthogonal, it is straightforward to show that Tr M µM T ν = δ µν (74) Suppose our tensor includes D irreps, all the js take D different values. Then there are D different on-leg singlets, and the matrix τ 0 is D 2 × D 2 ; thus, µ = 1, ..., D 2 and we have D 2M µ matrices. They act on the D dimensional space spanned by the D linearly independent on-leg singlets |0 (j) . These matrices form a D 2 linear space; we have shown thatM µ is an orthonormal set of D 2 matrices within this space, and thus it is an orthonormal basis and theM µ span the whole space of D × D real matrices.
The row transfer matrix and the norm thus take the formŝ and ψ|ψ = Tr Ê N = Tr E. Spectrum of the flux-free transfer matrix We have used the fact that each leg ofτ 0 forms a singlet |0 (j) ; however, recall the symmetry properties of the tensor A out of which the transfer operators were constructed, and the Gauss law satisfied by its four legs (46): j r ⊗ j u ∼ j l ⊗ j d . This implies that (j l ⊗ j r ) ⊗ (j d ⊗ j u ) must contain the single representation: the horizontal representations and the vertical ones must be such that can fuse to a singlet together. As a consequence of that, elements of (τ 0 ) j l ,jr;j d ,ju whose indices do not satisfy it must vanish. This splits the matrix (τ 0 ) j l ,jr;j d ,ju into separate blocks which can be separately diagonalized, implying similar block structure of the V matrices as well, splitting theM µ operators defined in (73) into different sets.
First, consider the so-called zero blockB 0 in which j l = j r as well as j d = j u . The elements of this block will be linear combinations of products of horizontal and vertical on-leg singlet projectors, TheM µ operators derived from this block will be diagonal in the space of singlets; the block (τ 0 ) jj;j j is a simple symmetric matrix, diagonalizable by the orthogonal block V jµ , using which we obtain the diagonal operatorŝ The next blocks are responsible toM µ which are offdiagonal in the singlet space. In the U (1) case, for example, we will have blocks for which j l − j r = j u − j d = ±k (for any integer k allowed by our tensors) Let us choose, in our U (1) example, to include one copy of each irrep |j| ≤ J (J may also be infinite). The matrix (τ 0 ) j l ,jr;ju,j d will have dimension of (2J + 1) 2 . The zeroth block ofτ , B 0 = (τ 0 ) j,jk;j ,j , will be a 2J + 1 dimensional matrix (since there are 2J + 1 possible on-leg singlet states). The blocks B k = (τ 0 ) j,j∓k;j ,j ±k will each be 2J + 1 − |k|) dimensional (counting the number of j values allowing for j −k and j +k values which agree with |j| ≤ J), from k = ±1 until k = ±2J -altogether 2J + 1 blocks whose dimensions add up, properly, to the right nally, since (τ 0 ) j l ,jr;j d ,ju is a symmetric matrix, we obtain that B k = B T −k , and write down the matrix in the block form (where the headers of the rows and columns denote the type of operators they connect with). This matrix can be easily blockwise diagonalized, involving the diagonalization of J +1 different blocks. Similar forms can be written also for other gauge groups (later on, we will work out a detailed example for the Z 2 case).
Before moving on to the contraction of Wilson loops, we shall consider some simple illustrative cases of norm computation, regardless of the gauge group. First, assume that all the blocks but the zeroth one vanish, and, on top of that, that the zeroth block is diagonal, that iŝ -all the relevantM µ operators are projectors (the other ones do not contribute since they are associated with zero eigenvalues). Then, it is easy to see that the transfer matrix iŝ E = j λ N j |0 (j) 0 (j)|⊗|0 (j) 0 (j)|⊗· · ·⊗|0 (j) 0 (j)| (82) Then, the eigenvectors are product vectors of the same representation, w j | = 0 (j)| ⊗ · · · ⊗ 0 (j)| with eigenvalues ρ j = λ N j , and the norm is Next, we keep the off-diagonal terms of the zeroth block zero, but allow for very small nonzero elements in the other blocks -that is, significantly smaller (in absolute value) than the diagonal terms of the zeroth block. If D irreps participate in our state, we havê M µ = |0 (j µ ) 0 (j µ )| for µ = 1, ..., D, with eigenvalues |λ 1 | ≥ ... ≥ |λ D | > 0; while for some K > D, |λ µ+1 | ≥ ... ≥ |λ µ+K | > 0 and there is some 1 ≤ L ≤ D for which |λ µ+1 | |λ L |. Then one may use perturbation theory to find the spectrum ofÊ. The nonperturbed part is · · · ⊗ |0 (j µ ) 0 (j µ )| giving rise to zeroth order eigenvectors as before, with corrections which are product vectors as well. Now allow for nonzero weak off diagonal elements in the zeroth block. Perturbation theory is still valid, keeping our eigenvectors close to product states along the row.
In fact, as long as the diagonal terms of the zeroth block are significantly stronger (in absolute value) than the rest of the τ 0 elements, this argument holds. As these other terms get larger and larger, the perturbative description loses its validity and the eigenvectors get farther from being product states along the row.
This may be interpreted as the lack or the presence of long-range order: the farther we are from product states along the row, the longer ranged order we have. Since confinement has to do with disorder [45], we find here the first hint to detecting area law from the transfer operators. As we shall see later on, indeed, the closer the transfer matrix eigenvectors are to product states, the closer we are to an area law of the Wilson loop.

V. A TALE OF TILING: CONTRACTING WILSON LOOPS
After having computed the norms, we move further to the contraction of Wilson Loop expectation values, which first requires studying further local ingredients: the fluxcarrying transfer operators.

A. Flux carrying transfer operators
Consider the transfer operators associated with sites carrying a straight flux line -that is, a group element operator U j (or U j † ) acting on either the horizontal or vertical direction, computed using (60) Using the symmetry conditions (44) as well as the transformation properties of the group element operators (31), we obtain that for every g ∈ G, θ l g ⊗θ †l g (85) (see Fig. 7(a)). That is, T → J M N maps from a total JM | on both ingoing legs (with respect to θ l g ⊗θ †l g ⊗ θ d g ⊗θ †d g ) onto JN | with respect to θ r g ⊗θ †r g on the outgoing horizontal leg and a singlet with respect to θ u g ⊗θ †u g on the outgoing vertical leg (see Fig. 8(a)). As in the flux-free case, that will have implications on the structure of the T → J M N operators.
Furthermore, the transfer operators T → J M N form a multiplet for each J, whose elements are mixed by the transformations. There is no problem with that, because in the contraction of the Wilson loop we sum over the M, N indices (matrix product and tracing of the U matrices). As usual, in the Abelian case the multiplets are trivial and contain one operator only, allowing us to give an intuitive illustration. For example, let us consider U (1) with the fundamental representation j = 1; there, the transformations take the simple form For the inverse horizontal flux line, one obtains (see Fig. 7(b)) -the difference from the right going flux is not very big, and has to do mainly on the opposite flux orientation: g instead of g −1 appears in the transformation, and the beginning index M is now associated with the right side rather than the left (similarly, N with the left rather than the right), since the flux goes backwards. This corresponds to transposition, and since the representations are unitary, D j nm (g) = D j mn (g −1 ) -i.e., the conjugate representation J. As a result, we denote the input of both legs as JN and the output of the right leg as JMvectors with a conjugate transformation rule (see Fig. 8(b)).
In the vertical direction, we have -the ingoing legs form a combined singlet, while both the outgoing legs, separately, belong to the J representation (one regular, one conjugate) -see Fig. 9.

B. Tiling the loop and projecting onto smaller spaces
Do we need to use all the elements of the transfer operators for the Wilson loop contraction? The answer is no; we can ignore some of them in the computation, while tiling the different building blocks together, thanks to the local symmetry and the special properties it enforces on the states and the transfer operators, just like we did in the case of the norm. As discussed, each of the local transfer operators used for the contraction, either with or without flux, can be seen as a map between the two ingoing legs to the two outgoing ones. While the ingoing legs form together a multiplet vector of the group, the output is a product of two separate multiplet vectors on the two outgoing legs (see Fig. 6(b), 8 and 9(b)). The numerator of the Wilson loop expectation value requires a particular tiling of the transfer operators, closing the loop. Since the output to each direction forms a multiplet vector, this will also be the input of the neighbouring transfer operators in the outgoing directions, and we can restrict all our transfer operators by cutting off all the input options that could not be realized within the Wilson loop tiling. This is done in a very similar way to what did in the norm computation, where we definedτ 0 (67) instead ofT .
Since our system is translationally invariant, let us identify the lower left corner of the loop with the lower left corner of our system. Let us consider the numerator of the expectation value of the Wilson loop: Where the trace is on both directions, assuming periodic boundaries (similar results may be easily derived for open boundary conditions); the J, M, N indices of the flux-carrying transfer operators have been omitted for simplicity, but it is assumed that they all carry the same irrep J (otherwise it would make no physical sense) and that the M, N indices are properly connected and summed over along the loop. Using the mapping properties summarized in figures 6(b), 8 and 9(b), we can write on each of the outgoing legs its output representa-tion -0, J or J for the conjugate representations used in the backwards fluxes cases. This immediately determines onto which inputs the transfer operators should be projected. Note that the lower right and both upper corners do not seem right in the equation above; nevertheless these are the right ingredients to be used, as explained below.
The tiling is composed of the following ingredients: Just like in the case ofτ 0 compared withT , these newly introduced operators contain less tensor elements and simplify the contraction of the Wilson loop, C. The decay of Wilson loops: is an area law possible?
Now we have all the ingredients required for the computation of a Wilson loop whose dimensions are R 1 × R 2 , and compute it using row transfer matrices, by contracting first in the horizontal direction, within an N × N system with periodic boundary conditions (torus).
We denote the transfer matrix corresponding to the first row we contract (the one containing the lower edge of the loop) by It is very similar to the MPS expression used for computing correlation functions (55) with one major difference. Due to the local symmetry, in between the two rows closing the loop, we need to use a different transfer matrix, E : the long range decay properties depend now two different transfer matrices, instead of one.
As stated in the beginning of this subsection, we have omitted the M, N indices and we assume implicit summation over them when contracting the loop. The Wilson loop contraction consists of the contraction of 2 (R 1 + R 2 ) indices, each taking dim (J) values -naively speaking, we would have to consider dim 2(R1+R2) (J) different contractions; however, the singular values are independent of these indices and depend only on the irrep J. Thanks to this symmetry, all the dim 2(R1+R2) (J) are equal, so it is enough to make one choice of the indices and multiply the result by dim 2(R1+R2) (J). This will be a perimeter-law term, however in the presence of an area law term it will not contribute in the large loop limit. Hence we focus below on computing for one particular choice of the indices.
Consider the diagonalization of the two transfer matrices which matter for the long range properties, Once again, we sort the eigenvalues in decreasing order, but in this case we do not care if the highest one is degenerate (but assume the existence of a spectral gap): for some integers K, K ≥ 1, Let us use this to compute the expectation value of the Wilson loop (107) in the thermodynamic limit N R 2 : (109) (we assumed that ρ 1 = ... = ρ K ; the generalization for the case of different phases is straightforward).
We further assume that the loop is large, that is -R 1 , R 2 1, allowing us to perform a similar simplification for E , and obtain that in the thermodynamic limit, for large loops, (This holds only if if this condition is not fulfilled, the vectors v j and w i | should not be seen as those corresponding to the highest eigenvalues, but rather as those with the highest eigenvalues for which this condition is satisfied. We assumed here that ρ 1 = ... = ρ K ; the generalization for the case of different phases is straightforward).
Assuming rotational invariance, we could repeat the same procedure by contracting the columns first, to obtain Both expressions must be equal; therefore, we deduce that But the more interesting question is whether ∂ρ 1 (R) /∂R = 0 or not. If the largest eigenvalue ofÊ (R) does not depend on R, we obtain that with some constantC: perimeter law decay of the Wilson loop (unless ρ 1 = ρ 1 ). On the other hand, an area law is possible if with κ > 0. Let us plug this expression into (111) and (112). We will obtain the equation for some constant C, and we and obtain, finally, for large Wilson loops, that if ρ 1 (R) ∼ Γe −κR , (117) -exactly the same form of (39), with W 0 = C KΓρ1 , κ A = κ and κ P = log ρ1 Γdim 2 (J) − κ. Therefore, we conclude that a perimeter law will be obtained if the largest relevant (in terms of accessible throughÊ b andÊ t ) eigenvalue ofÊ (R) is independent of R; an area law is possible if it depends on R exponentially. Why only possible? To see why this condition is necessary but not sufficient for the area law to hold, let us consider the following scenario.
Previously, we made the assumption that the eigenvectors of the flux-free transfer matrix should be close to product vectors in order to make an area law possible. We also know that the expectation value of the Wilson loop depends on the zeroth flux transfer operatorsτ 0 inside and outside the loop, and some other, flux-carrying transfer operators along the loop. Let us assume that we are, indeed, in a scenario in which the eigenvalues of the transfer matrix are close to product states. Denote as usual the highest eigenvalue of the transfer operator by λ 1 . Then the norm, for a large system, will roughly scale as λ N 2 1 : each site contributes a single power of λ 1 . This is the denominator of the expectation value formula. In the numerator, we will have a contribution of λ 1 for each site outside the loop; within the loop, it depends.
If the flux carrying transfer operators along the loop take us from the singlet subspace corresponding to λ 1 to that of another eigenvalue -denote it by λ -we will have a contribution of λ for each of the sites within the loop, and the Wilson loop's expectation value will scale as (λ /λ 1 ) A where A is the area of the loop (E (R) ∝ λ R ). However, if the flux carrying transfer operators do not take us to another singlet subspace with a different eigenvalue, we will not have an area dependent contribution. In this case, the largest eigenvalue ofÊ (R) depends exponentially on R (through λ R ) but an area law is not obtained, which shows us why this condition is necessary but not sufficient.
On the other hand, if the eigenvectors ofÊ are far from product vectors, which means they are governed by some collective, long range effect, we cannot have areadependent contributions at all.

VI. ILLUSTRATION: THE Z2 CASE
To conclude and illustrate our discussion, we will show an explicit example, where the gauge group is Z 2 . In this case, the group Hilbert space on each link is two dimensional, with representations labelled by j = +, −, which can be simply seen as spins. The group element operators are Hermitian, U = U † = X, and invert the spin, and the group operations Θ (no difference between left and right in Abelian groups) are the identity operator as well as Gauge transformations are given bŷ (120) We would like to consider the most general PEPS with translational and rotational invariance, with physical spaces containing all the irreps and virtual ones containing a single copy of each irrep (minimal construction -as explained above, to consider real physical scenarios one will most likely have to generalize in a straight forward manner and add more copies, as was necessary in the Z 3 demonstration of Ref. [42]). Thus, the physical and virtual spaces will be the same, two dimensional spin-like spaces spanned by the representation states |± . The state will be parametrized by the tensors A st lrdu , with s, t, l, r, d, u = ±. The most general construction satisfying these conditions is given by and the rest of the elements, which violate the symmetry, vanish. If we consider the |+ states as flux free states, and the |− as flux carrying, we can interpret α as the amplitude of having no fluxes going through the site, β as the amplitude of corner flux, γ -of straight line fluxes and δ -two intersecting flux lines.
Here we will be interested in the properties of the transfer operators constructed for such states, and the computation of the Wilson loop expectation value.

A. The transfer operators
The transfer operatorT may be simply built using (58) and (59).
Let us identify the elements of the vector space spanned by the double legs of the transfer matrix. The on-leg transformations here admit the simple form θ ⊗ θ † = Z ⊗ Z for the only group element which is not the identity; Since there are two irreps, we will have two onleg singlets, as well as two non-singlets, Where the new notation introduced in the two equations above factorizes the on-leg Hilbert space into the product of two spin spaces; one detects whether the state is an on-leg singlet (s) or not (n) and the other labels the two states within each of these options by ↑ and ↓.
Using these states, we can write down all the relevant transfer operators and their reductions. For example, where the block structure is clearly seen; the first one is the zeroth block, mixing only projection operators. It depends on α, γ, δ -the amplitudes for which fluxes do not change directions, and thus the representations are not changed horizontally and vertically on the state, and the on-leg singlets are not flipped on the transfer operators. The second block, where the representation / singlet change, depends on β -the turning (corner) flux amplitude. Furthermore, as the parameter γ has to do with straight flux lines going through the site, we expect that the larger it gets, the farther theM µ operators derived from the zeroth block are from projection operators, and the farther we are from an area law; indeed, as we see, it appears on the off-diagonal terms of the zeroth block, and when γ = 0 theM µ operators of the zeroth blocks are projectors.
This matrix can be easily diagonalized as in (71), with the eigenvalues (not necessarily in descending order -this depends on the values of the parameters): with the diagonalizing matrix |↑ ↑| ⊗ |s s| u 11 (α, γ, δ) u 12 (α, γ, δ) 0 0 |↓ ↓| ⊗ |s s| u 21 (α, γ, δ) u 22 (α, γ, δ) 0 0 |↑ ↓| ⊗ |s s| 0 0 Using all that, we obtain the operatorsM µ as defined in (73), where Π ↑ = |↑ ↑| and Π ↓ = |↓ ↓|. Since V is orthogonal, they form an orthonormal basis as in (74).M 4 is irrelevant, since λ 4 = 0; the |s s| is also irrelevant since it multiplies everything, and hence we will omit it and refer to the operatorsM µ as two dimensional. Note that as expected the first two ones,M 1,2 , having to do with the zeroth block, are diagonal, while the other ones are not. Similarly, we can compute and write down the other relevant matrices. Note that since the fluxes have no orientation in our case,τ → =τ ← ≡τ − andτ ↑ =τ ↓ ≡τ | . We thus require only six rather than eight further matrices. The first iŝ connecting operators acting on the non-singlet subspace in the horizontal direction with ones acting on the singlet space in the vertical one. The same block structure is apparent; the first block is a generalization of the zeroth block -still only connecting projection operators, though acting on different spaces, and the second block changes the representations. As in the τ 0 case, the parameter γ is the one "spoiling" the area law: all the amplitudes ofL µ operators which do not change the on-leg singlet eigenvalue subspace are proportional to it. One it is set to zero, when crossing a flux line the subspace will change. We can formally perform a horizontal-vertical singular value decomposition and obtain an expression of the form τ − = µ η µKµ ⊗L µ . Since the horizontal operators act only within the non-singlet subspace and the vertical ones only within the singlet subspace, we can representK µ andL µ by two dimensional matrices.
τ | is simply obtained by transposition, Finally, let us consider the transfer operators of the four corners. We begin with the lower left corner τ = lr / du Π ↑ ⊗ |s n| Π ↓ ⊗ |s n| σ + ⊗ |s n| σ − ⊗ |s n|  where in both dimensions we get a singlet input and obtain a non-singlet output. Here, after performing the singular value decomposition, we will also use two dimensional operators acting only on the "spin space" since this corner operator connects to the right s/n subspaces. The other corner operators arê Note that all the elements of the corner operators are proportional to either β or β, which is expected since β is the corner parameter, and it would be impossible to close a loop in its absence.

B. Analytical example
Let us set, for simplicity, γ = 0. Consider τ 0 (125) and theM µ operators derived from it (128). Let us set γ = 0; then we simply haveM as well as The choice of γ = 0 sets all the zeroth blockM µ operators to projectors onto orthogonal states, and the flux-free transfer matrix from (75) takes the form If we further assume that |β| |α| , |δ| we find ourselves in the perturbative case discussed above, and may use perturbation theory for finding the eigenvectors ofÊ. The zeroth, unperturbed part is in the first row of (137), from which we find two approximate, zeroth order eigenvectors, with zeroth order eigenvalues -which are the two highest ones. Let us set, without losing generality, |α| > |δ| (one can easily invert that in the following discussion). The leading order corrections to the eigenvalues will be second order (∝ |β| 8 ) and to the eigenvectors will be of the first order (∝ |β| 4 ); we shall neglect them both. The norm of the state is then That is, the torus is tiled with N 2 sites, each contributing a factor of |α| 2 to the norm. Let us now move on to the flux carrying transfer matrices. Looking at the straight flux onesτ − (129) andτ | (130), we see that our choice of γ = 0 sets the zeroth block to zero. This implies that they will flip the local incoming spins in both directions -in particular in the direction orthogonal to the flux; i.e., the eigenspace ofτ 0 out of the loop will be connected to the orthogonal one within the loop, eventually to give rise to an area law, unless |α| = |δ|. We see thatτ (ignoring the n, s space for the reasons explained above) -inverting the spins in the orthogonal direction to the flux lines, that is, changing indeed from the α to the δ sector and vice versa. For the corners we get Let us consider the action of the lower row of the Wilson loop,Ê b (R 1 ) on the input state w 1 | with the highest eigenvalue, identifying without loss of generality, as usual, the origin of the torus with the lower left corner of the loop. We get (143) where the omitted terms either annihilate w 1 | or are of negligible magnitude.
The leading terms of the output vector w 1 |Ê b (R) are product vectors, with ↓| entering the loop and ↑| out of it.
The two spins which are on the loop's boundaries are either flipped or not, depending on the particular configuration from the list above. We get for an even R 1 (144) and for an odd R 1 We move on to the intermediate rows, witĥ where, once again, the terms not included are either small enough or annihilate the input vector. The highest eigenvalue (in absolute value) is -exponential in the distance R 1 , just as speculated in (114), with Γ = |α| 2(N −1) β 2 δ 2 , and string tension κ = −2 log δ α -predicting an area law behaviour. This eigenvalue is four-fold degenerate (in absolute value). Denoting by |x = ±1 the eigenvectors of σ x , with eigenvalues ±1, we get the four eigenvectors, (148) Note that since the transfer matricesÊ andÊ (R 1 ) are hermitian, |v i = |w i and |v i (R) = |w i (R) .
Implying that for an even R 1 Giving rise to, for an even R 1 and for an odd one We are finally ready to obtain the Wilson loop expectation value using the procedure of section V C. We will have to slightly modify it, since in our case the highest eigenvalue ofÊ is only degenerate in absolute value; for large loops in the thermodynamic limit we thus modify Eq. (110) to where p = even,odd is the parity of R 1 . One can already clearly see the area and perimeter dependent parts. The only thing left to do is to complete the computation of the sum, where four different cases have to be considered, corresponding to the parities of R 1 , R 2 . It is straightforward to see that if the area is even (three of the four cases), the resulting number is 8, while if the area is odd, the result is 8Re αβ |αβ| 4 . Altogether we obtain, for large loops in the thermodynamic limit, for the |β| |δ| < |α| and γ = 0 case, that The Creutz parameter (38) is nothing but the string tension, We see that we have an area law, or a confining phase, as long as |δ| = |α|. While we excluded an equality in our arguments above, indeed we will have no area law if these two parameters are equal: then, the eigenvectors ofÊ between which the fluxes transfer will have the same eigenvalue which does not allow for an area law, in full accordance with our general discussion. If we switch γ on, it will have two effects: one will contaminate the eigenvectors of the transfer matrixÊ, taking them farther from product vectors until the area law is broken, as well introduce terms in the flux-carrying transfer matrices that do not change the eigenvalue sector ofτ 0violating another area law criterion.

C. Numerical examples
We will now present a few more examples which are computed numerically, using exact contraction, on a torus with size N 1 = 8 × N 2 = 100. We considered different choices of parameters to demonstrate different behaviours; for each, we computed expectation value of the Wilson loop for several large loops. We extracted the parameters κ A and κ P as follows: using the expression (39) for a Wilson loop, we may define a function of R 2 depending on R 1 as a parameter, f (R 2 ) = − log W (R 1 , R 2 ) = f 1 (R 1 ) R 2 + f 0 (R 1 ) (159) It is a linear function, which intersects with the vertical axis at f 0 (R 1 ) = κ P R 1 − log W 0 (160) whose slope is In the case of a perimeter law, the slope function will be constant, f 1 (R 1 ) = κ P and when plotting f (R 2 ) for different R 1 values, parallel lines will be obtained. In the case of an area law, the lines will have different slopes. Thus, κ A and κ P may be extracted by performing linear fits to the functions f 1,2 (R 1 ). Moreover, we have extracted the Creutz parameter too. The first set of parameters we examine is α = 1, β = 0.1, γ = 0, δ = 0.95. This choice is within the perturbative class studied above. It shows an area law, as can be seen from Fig. 10 and the Creutz parameter χ = κ = −2 log δ α ≈ 0.1025 (as shown in Fig. 11). The expected exponential dependence of the eigenvalues ofÊ (R) is demonstrated in Fig. 12.
Next, let us consider another example which lies within the perturbative regime: α = 1, β = 0.1, γ = 0, δ = 1. Here still γ = 0 and β is very small, so the eigenvectors ofÊ would be product vectors, hence satisfying the first criterion for an area law. However, the eigenvalues ofτ 0 are degenerate, implying no area law (the second criterion is violated). The perimeter law is clearly shown in Fig.  13, and, as as one can see in Fig. 14, the eigenvalues of E (R) have no dependence on R.
Finally, we consider a completely different case, where α = 0.1, β = 0.1, γ = 1, δ = 0.3. For this choice of parameters, the previous perturbative treatment is not valid. The eigenvalues associated withτ 0 are (163) -here, too, the most significant contributions are from the zeroth block with diagonal operators (the first two); however, they are far away from being projectors, hence we do not expect the eigenvectors ofÊ to be anywhere close to product vectors.
Let as also consider the straight flux carrying transfer Figure 10. The α = 1, β = 0.1, γ = 0, δ = 0.95, which lies within the perturbative class discussed above, clearly shows an area law. It can be seen qualitatively on the top, where − log W (R1, R2) is plotted as a function of R2 for three different values of R1 -resulting in three non-parallel lines. And if it is hard to detect the different slopes on the top, the middle figure shows it more quantitatively: the slope function f1 (R1) ≈ 0.126R1 + 9.1078 has a nonzero slope κA ≈ 9.1078, and its intersection with the vertical axis is κP ≈ 9.1078, the slope of the function plotted on the bottom, f0 (R1).
operatorτ − and toτ | . We find the singular values η 1 ≈ 0.4472, η 2 = 0.02, η 3 = η 4 = 0, (164) Figure 11. Computation of the Creutz Parameter χ (R1, R2) for α = 1, β = 0.1, γ = 0, δ = 0.95, for different values of R1 and R2. As can be seen, the values converge to the predicted value (thanks to the validity of the perturbative treatment in this parameter regime) of −2 log δ α ≈ 0.1025. Figure 12. In the perturbative case worked out analytically, α = 1, β = 0.1, γ = 0, δ = 0.95, the highest eigenvalue of the intermediate transfer matrixÊ (R) depend exponentially on the width R, as can be seen from the logarithmic plot given above, where the two highest eigenvalues (in absolute value, both degenerate in this case) are plotted for all values of R. The symmetric shape is due to the finiteness of the system (N = 8 in this case). For R ≤ 4, the eigenvectors corresponding to highest eigenvalue connects with the right input state, while for R ≥ 4 the next ones are relevant -all due to the symmetry. Also shown is a linear fit, computed with respect to the parameters predicted using the perturbative treatment.
associated, in the flux direction, with the operatorŝ and, in the direction orthogonal to the flux, with the operatorŝ which imply that even if our eigenvectors were product Figure 13. The α = 1, β = 0.1, γ = 0, δ = 1 case does not allow for an area law because of the degeneracy in the eigenvalues of τ0 corresponding to projection operators. The state shows a perimeter law, which can be seen qualitatively on the left, where − log W (R1, R2) is plotted as a function of R2 for three different values of R1 -resulting in three parallel lines. Quantitatively we see in the middle, where the three slopes of the three lines are plotted, that they are equal: f1 (R1) = κP ≈ 9.2103 is a constant function (κA = 0). On the right we see the fit of f0 (R1) ≈ 9.2103 (R1 − 1). Figure 14. The α = 1, β = 0.1, γ = 0, δ = 1 case does not allow for an area law because of the degeneracy in the eigenvalues ofτ0 corresponding to projection operators. This is also manifested by the fact that the eigenvalues of the intermediate transfer matrix,Ê (R), are completely independent of the distance R, as illustrated here.
vectors (which they are not), the most prominent contribution, coming from η 1 , would be diagonal in the subsector (as seen fromL 1 ). Therefore all our area law criteria are violated. Indeed, this set of parameters show a perimeter law decay of the Wilson loop, as can be seen in Fig. 15, in the zero Creutz parameter (see Fig. 16) and in the eigenvalues ofÊ (R) which are independent of R (as shown in Fig. 17).

VII. SUMMARY
In this work we have seen how local properties of two dimensional lattice gauge theory PEPS, manifested in their transfer operators (on-site) and matrices (rows) simplify their contraction and dictate their long-range, Wilson loop behaviour. We have related the area law with transfer matrices whose eigenvectors are product vectors -that is, a product of local contributions of the transfer operators on each side, manifesting the lack of long-range order, as expected for a disordered, confining phase. The perimeter law, appearing in ordered phases, has to do with non-product eigenvectors, where the separate sites contribute in a correlated, long-ordered man-  ner. These results may be used for detecting phases of PEPS used for pure gauge theory studies, and for the design of PEPS used as ansatz states for such scenarios.
One possible extension is the inclusion of dynamical matter -which is different from the current work both in the mathematical sense (different structure of the tensors, implying different symmetry properties) and the physical one (in that case, at least with fermionic matter as in conventional standard model scenarios, the Wilson loop does not serve as an order parameter for confinement any more). This could possibly connected with the formalism of gauged Gaussian fermionic PEPS [31,32] which can be contracted using sign-problem free Monte-Carlo techniques [33,37] both for the study of further examples and application to physical models of interest.
Another important and relevant generalization is the extension to higher dimensions, where further geometry arguments have to be taken into account, potentially containing many further interesting physical and mathematical properties.