Variational Monte Carlo simulation with tensor networks of a pure $\mathbb{Z}_3$ gauge theory in (2+1)d

Variational minimization of tensor network states enables the exploration of low energy states of lattice gauge theories. However, the exact numerical evaluation of high-dimensional tensor network states remains challenging. By combining gauged Gaussian projected entangled pair states with a variational Monte Carlo procedure, we are able to efficiently compute ground state energies. We demonstrate the method for a pure gauge Kogut-Susskind Hamiltonian with a $\mathbb{Z}_3$ gauge field in two spatial dimensions. The method provides an inherent way to increase the number of variational parameters and can be readily extended to systems with physical fermions.


I. INTRODUCTION
Tensor network states, especially matrix product states (MPS), have changed our understanding of solid state systems dramatically. Describing states with an arealaw entanglement, i.e. ground states of local, gapped Hamiltonians [1,2], MPS provide an ansatz class for a wide range of problems due to their favorable numerical scaling. Instead of an exponential scaling, MPS algorithms scale polynomially with the system size. The computational power in combination with a solid analytical understanding allowed a variety of applications, including ground state searches [3,4] and the description of dynamics of many-body systems. Similar studies have been performed with tensor networks in two spatial dimensions, projected entangled pair states (PEPS) [5].
Motivated by the success of tensor networks in condensed matter physics, such methods have been generalized and applied to particle physics problems too, in particular to lattice gauge theories (LGTs) [6]. Gauge theories appear in many fundamental physical contexts, e.g. the standard model of particle physics, where gauge fields act as force carriers. In particular, it includes Quantum Chromodynamics (QCD), the theory of the strong nuclear force, which, as a non-Abelian gauge theory [7] has a running coupling. In QCD, asymptotic freedom [8] gives rise to asymptotically weak couplings for high energy scales (e.g. collider experiments), and therefore perturbation theory could be used in these physical regimes. On the other hand, low energy QCD is a strongly coupled model, requiring non-perturbative treatment.
One approach to regimes where non-perturbative methods break down are lattice gauge theories. They provide a gauge invariant regularization of gauge theories, discretizing either spacetime [9] or only space (leaving time continuous) [10]. Simulations based on hybrid Monte Carlo [11,12] have given many interesting insights into the physics in the non-perturbative regime. While having been extremely successful and fruitful for static studies (such as studies of the hadronic spectrum), this method faces two major difficulties. First, the inability to directly observe time dependent phenomena in Wickrotated, Euclidean spacetimes, as done in this context; the second is the well known sign problem [13] which appears in scenarios with finite fermionic chemical potential, where the statistical interpretation allowing one to perform Monte-Carlo sampling breaks down, blocking the way to important phases of the QCD phase diagram [14].

In (1+1)d, MPS have been very successful describing
LGTs (see Ref. [6] and references therein). In higher dimensions, MPS are generalized to PEPS, whose contraction is in general very costly. This hinders the application of variational PEPS algorithms in higher dimensions, although a first numerical study for a pure gauge theory has been recently presented in [15]. Earlier numerical studies used less general tensor networks for two dimensional lattice gauge theories, either purely gauge [16] or including fermions [17]. In contrast, analytical approaches have developed faster, with the formulation of gauge invariant pure gauge PEPS [18], and more general gauging mechanisms including matter for arbitrarily dimensional PEPS [19,20].
In these works, the global symmetry of a matter-only PEPS is lifted to a local one by introducing a gauge field, in a way analogous to minimal coupling. The latter gauging method has been used for the construction of gauged Gaussian fermionic PEPS [21,22], where the matter state to be gauged is a free (Gaussian) fermionic state, in a manner analogous to minimal coupling of a Hamiltonian [23]. The restriction to this subclass of PEPS enables the efficient contraction of the states with Monte Carlo techniques [24]. Since the sampling probability of the algorithm depends only on the norm of the state, the Monte Carlo algorithm cannot suffer from the sign-problem. Furthermore, the construction allows for a natural and efficient extension to higher bond dimensions which is numerically very expensive in general PEPS calculations. However, until now, these states have only been used to compute observables of a toy model without an energy minimization. It is a priori unclear whether the energy minimization of a real Hamiltonian will converge to the ground state energy and whether the expressability of the Ansatz is sufficient to capture the relevant physics of a theory.
In this paper, we present the application of gauged Gaussian PEPS [20][21][22]24] as ansatz states in a variational Monte Carlo (VMC) procedure [25][26][27]. We apply the algorithm to a Hamiltonian pure Z 3 gauge theory [28] and make explicit use of the possibility to extend the Ansatz efficiently by adding more layers of virtual parameters.
The Z 3 theory is a relatively simple (2+1)d theory, but it is known to exhibit a (first-order) phase transition between a confining and non-confining phase, and thus constitutes a non-trivial testbench for the ansatz [29]. Furthermore, extensive Monte Carlo studies have been performed on Z N theories, which allow us to benchmark our results against known results [30]. Our goal is to demonstrate the expressibility of the ansatz presented in Ref. [24] and how it can be applied to study gauge theories. Adding more layers to the construction is essential to improve convergence, especially in the low coupling regime of the theory. However, precisely locating the phase transition remains challenging, even with an increased number of layers. The main obstacle is the expensive evaluation of a Pfaffian that appears in the calculation of the electric energy. Thus, it has to be calculated in every Monte Carlo step during the energy minimization.
The rest of the manuscript is structured as follows: In sections II and III, we introduce Z N gauge theories and construct our ansatz states. These states are minimized with the numerical methods described in section IV. The numerical results are presented in section V. Finally, we conclude in section VI.

II. HILBERT SPACE OF ABELIAN LATTICE GAUGE THEORIES
In a Hamiltonian lattice gauge theory, space is discretized and represented on a lattice while time remains continuous [10]. This is in contrast to the action formulation, where both space and time are discretized [9]. The (fermionic) matter of the theory resides on the vertices x of a lattice, and the interactions are mediated by gauge fields, whose quantum Hilbert spaces reside on the links (cmp. Figure 1). In the following, we will focus on Abelian lattice gauge theories with finite gauge groups (Z N ), without dynamical matter, i.e. pure gauge theories. We will consider a two-dimensional L × L lattice with periodic boundary conditions. Thus, the only degrees of freedom of the theory reside on the links.
One problem of numerically simulating a lattice gauge theories with compact Lie groups (even the Abelian U (1)) is the infinite dimension of the Hilbert spaces on the links. This can be approached by truncating the local Hilbert spaces, either by introducing a cutoff to the electric field, allowing one to restore the full theory  by extending the cutoff [21] or integrating over an extra dimension [31][32][33], or by sampling group elements [28] from the gauge group, which form a subgroup. Due to the construction of our states (see section III), we chose the second approach, i.e. instead of simulating the full U (1) theory, we consider a Z N subgroup that serves as an approximation for U (1). As described in Ref. [28], the N → ∞ limit of Z N reproduces U (1), and hence Z N lattice gauge theories flow, in the large N limit, to compact QED [34], a lattice gauge theory with U (1) symmetry.
We write the Hamiltonian of a pure Z N gauge theory as where = (x, i) is a link on the lattice emanating from vertex x horizontally(i =ê 1 ) or vertically (i =ê 2 ) and p is a plaquette [28]. The indices p j refer to one of the four links of one plaquette as indicated in Figure 2. The terms H E and H B are referred to as electric and magnetic part of the Hamiltonian, respectively [10].
The operators in (1) obey the Z N algebra given by Operators that act on different links commute with each other. The Hamiltonian (1) is invariant under the action of the local unitary operators The links are addressed according to their vertex x and their direction right (r) or up (u). This local gauge invariance implies that Θ(x) commutes with the Hamiltonian on each site Due to the generators of local symmetry (given in (3)), we know that the physical states of the system obey the symmetry Equation (5) holds since we do not consider static charges in this work. Given the Z N group, we define a set of group element states |q( ) labelled by integers q = 0, ..., N − 1, which span the local gauge field Hilbert space on link . They correspond to group elements with the discrete angles φ(l) = qδ (δ is defined in (2)). The group element states form an orthonormal basis for the local Hilbert space q|q = δ q,q .
These states are eigenstates of the Q operators, with They are lowered by the P operators, periodically:

III. PEPS CONSTRUCTION WITH ABELIAN SYMMETRY
Products of local group element states define the configuration of gauge fields on the lattice. Such product states, |G = ⊗ |q( ) form an orthonormal basis, using which we can expand every state in the gauge field Hilbert space: where the sum runs over all possible gauge field configurations on the links and Ψ(G) is a gauge field dependent wave function of the configuration G. This expression is a special case of the more general formulation presented in [24], where Ψ(G) can be a quantum state of the dynamical (fermionic) matter, |Ψ(G) , instead of the wave function we have in our current pure gauge case. Not every state that can be expressed with (8) is physically relevant, i.e. fulfills the local symmetry (3). Thus, the wave function Ψ(G) has to be chosen such that the full state |Ψ obeys the correct symmetries. Additionally, the state that we pick should allow for efficient numerical calculations of observables and gradients. Following the general construction in [24], we choose a gauged Gaussian projected entangled pair state (GGPEPS) as an ansatz. For details and further motivation, we refer to Refs. [21,23].

A. Construction with a single layer
Following the idea of a PEPS construction, we build the GGPEPS out of local constituents which help us to impose the symmetry. The local parts are entangled during the construction to form the final wave function.
The elementary building blocks for the wave function are auxiliary (or virtual) fermionic modes that are attached to each outgoing and ingoing leg of each vertex of the lattice. They are chosen to be fermionic to enable a consistent coupling to fermionic matter which obeys the correct statistics [24]. Although, for the description of a pure gauge theory, the coupling to matter is not necessary.
The construction of a GGPEPS consists of three essential parts (cmp. Figure 3). First, the fiducial operators A(x) create virtual fermionic states out of the modes associated with each site. They are constructed in a way that guarantees virtual gauge invariance (used in general PEPS constructions for imposing global symmetries). This step of the construction can be readily extended to include more virtual fermions, in a similar spirit that the bond dimension of a PEPS can be increased. The details of the construction with multiple layers are given below. Then, some of the virtual modes on each site are rotated with respect to the physical gauge fields of the theory, in a particular way that lifts the virtual symmetries to physical ones [24]. This is done by gauging operators U G acting on the virtual fermions and controlled by the gauge field configuration. Finally, the pairs of virtual fermionic modes on the two sides of each link are projected onto maximally entangled states by projection operators ω . That contracts the state from its local constituents and introduces correlations to the state.
The wave function can thus be written as where the products are over all links of the lattice and |Ω v is the fermionic Fock vacuum. In the following, we will treat the three main components of the construction A, U G , and ω in more detail, and see how to make sure that Ψ(G) obeys the right symmetry properties. Furthermore, aiming at an efficient computation of the wave function, we would like it to be Gaussian, and thus all its constituents will be Gaussian too.
On each vertex x of the two dimensional lattice, we define eight virtual fermionic modes, two associated to each leg -left, right, up and down. On each leg we label the two modes by ±, and sort them into two groups: a i = {l + , r − , u − , d + } (which we call the negative modes) and b i = {l − , r + , u + , d − } (positive modes). The modes obey the Dirac anti-commutation relation c(x), c † (y) = δ x,y and {c(x), c(y)} = c † (x), c † (y) = 0, where x, y are vertices on the lattice and c is a fermionic mode. We define the virtual electric fields with k ∈ {r, l, u, d} as well as the generator of the gauge transformation on the virtual degrees of freedom, This can be seen as a version of a Gauss law operator: the divergence of the virtual electric fields at the vertex. The staggering is introduced to accommodate the general case with physical fermions [24] (aiming at the problem of physical fermion doubling [35] which we do not encounter in the pure gauge case). It is taken care of already on the level of electric fields (cmp. (10)) and thus the rest of the equations can be stated without explicit reference to staggering. The fiducial operator A(x) which creates the modes out of the vacuum has to be Gaussian, and be invariant under transformation generated by G 0 (x). Hence, it is given by [21,23] where T ij is a 4 × 4 matrix containing all parameters of the ansatz. A is a Gaussian operator by construction, and one can easily inspect that since positive modes are only coupled to negative ones, the symmetry property is satisfied for every angle α, hence forming a U (1) parameterization. As such, it holds also for the Z N cases, with a discrete choice of angles. Due to other symmetry considerations (e.g. lattice rotation invariance), only two independent parameters in T ij of initially sixteen remain, y and z. They couple different modes in a given vertex: y couples right(up) and left(down) modes in a vertex, z couples modes that are building corners, e.g. right and up modes. The exact form of T and a motivation of the symmetries can be found in Appendix A.
For now, we will formulate the ansatz with eight virtual fermions per vertex. One set of eight virtual fermions is referred to as one layer. In a second step, we will enlarge the number of variational parameters by adding more layers, i.e more virtual fermions to the links. Each layer gets an independent set of parameters y and z. Increasing the number of layers is the analogue to increasing the virtual bond dimension in a non-fermionic PEPS.
In a second step, we entangle the virtual fermions on the links with physical gauge fields on the links. The gauging operator for a given gauge field configuration G takes the form where q( ) parameterizes the group element on the link in the configuration G. The local gauge transformation changes only the modes pointing up and right. Modifying the left and bottom modes as well would undo the gauge transformation due to the staggering. For a detailed overview of the gauging procedure in terms of PEPS operators, i.e. in graphical notation, we refer to Refs. [20,21,23].
In order to create more than a product state, we project the virtual, fermionic modes adjacent to each link onto maximally entangled states. The unnormalized pro- connect the left(upper) and right(lower) modes of neighboring sites. Here, Ω is the projector to the virtual vacuum on link andê i is the unit vector in direction i. Similar to the fiducial operators A, the projectors ω are Gaussian and commute among each other since they are products of fermionic modes on different links. The projectors link the virtual modes of one site with the virtual modes of the next site in horizontal and vertical direction, respectively. It is essential that the projectors are unnormalized since the norm of a state will serve as a transition probability between different gauge field configurations later.
Combining A, ω, and U G , we get the wave function in equation (9). Now, we can show that the construction is indeed gauge invariant and fulfills (5). We act with Θ(x) on |Ψ explicitly, on some given vertex x whereq are all gauge fields that are not affected by the gauge transformation, i.e. that are not adjacent to x. To shorten notation, we named the different links according to the labels defined in Figure 2. The third line is linked to the second one by a change of variables in q. The gauge invariance holds if Ψ(G) = Ψ(G ). We can write the wave function Ψ(G ) as where˜ are all links that are unaffected by the gauge transformation and Ω v is the vacuum of all virtual modes.
The notation of multiple signs shows the transformation for an even (top sign) and an odd (bottom sign) vertex at the same time. We used the invariance of the fiducial operator (13) at the last line. In order to transform the virtual electric field from the adjacent vertices x −ê 1 and x−ê 2 to vertex x, we use the invariance of the projectors ω ω x−ê1,1 e iδE0(x−ê1,r) = ω x−ê1,1 e −iδE0(x,l) All operators employed in the construction (A, ω, and U G ) are Gaussian operators. Since products of Gaussian operators are still Gaussian [36], the wave function Ψ(G) can be efficiently described with covariance matrices. As detailed in [24], there are multiple ways of combining the operators to covariance matrices. We choose to group the gauging operators and the projectors together into Γ in (G), a covariance matrix that depends on the gauge. The fiducial operators are summarized in a second covariance matrix D. The relation between the covariance matrices and the gauged ansatz state can be summarized as For further details about the formulation of Gaussian operators in terms of covariance matrices, see Appendix C.
The covariance matrices or parts of them allow the efficient calculation of the Monte Carlo transition probability (cmp. equation (25)).

B. Construction with multiple layers
Although the ansatz wave function with a single layer, i.e. two variational parameters, captures the high coupling regime very well, the low coupling regime is challenging for a single layer (cmp. Figure 5). Upon increasing the number of layers, the agreement between exact diagonalization data and the variational PEPS approach improves dramatically. In order to increase the number of variational parameters, we add more virtual fermions to the construction. Each layer carries an independent set of parameters, i.e. the matrix T in the fiducial operator A is different for each layer, while the states are coupled to the same gauge field. This ensures that all states fulfill the Gauss law. The virtual fermions of different layers on the links do not interact. The complexity of the computation scales linearly in the number of layers because the state can be contracted as independent layers of equally sized PEPS. Further details about the contraction and the changes to the calculation of observables are explained in Appendix B.

IV. COMPUTATIONAL EVALUATION
The ansatz defined above characterizes a family of states that depends on two parameters. In order to find the ground state of the Hamiltonian (1) for N = 3, we have to adapt the parameters such that the energy is minimized. By computing expectation values of observables and derivatives with respect to the parameters via sampling, we circumvent the unfavorable scaling of PEPS contractions. The variational Monte Carlo technique works in a two step procedure: first, the energy and the gradients are sampled for a given set of parameters α. In the second step, the parameters are changed α → α according to the gradients and a minimization algorithm.

A. Calculation of expectation values
The Hamiltonian (1) consists of two terms, the electric energy and the magnetic energy. Due to translational invariance of the states and the Hamiltonian, it is sufficient to calculate the energy of a single plaquette and a single link where n plaq = L 2 , n links = 2n plaq and L is the linear extent of the quadratic lattice (number of vertices). In the equation above, is a freely chosen link. If not stated otherwise, we choose the link at x = 0 in the horizontal direction. Calculating the magnetic energy is a special case of the expectation value of a Wilson loop. We define the Wilson loop operator as where C is an oriented, rectangular curve of length R 1 in the horizontal and R 2 in the vertical direction. The oper- ator Q is picked as is or daggered according to whether the link is traversed in the direction of the blue arrows (cmp. Figure 4) or against them. The Wilson loop operator does not only play a role for the calculation of the energy, but can be used as indicator for confinement in the theory (cmp. section V). Given the state defined in (8), the expectation value of a Wilson loop reads where the estimator F W (R1,R2) = ∈C exp(±iφ( )) is a complex number and the sampling probability is While the expression · is the expectation value of an operator, the expression · MC is a p(G)-weighted average over complex numbers. Since the norm of a state is always real and larger than zero, this formulation of a Monte Carlo procedure cannot suffer from the sign problem.
Using the covariance matrices defined in (20) in the formulation of Majorana fermions (cmp. Appendix C), we can write the squared norm of the wave function as It serves as the transition probability between different configuration states of the gauge field. In our Monte Carlo scheme, we use the Metropolis algorithm [37] with eq. (24) as a transition probability. In each step, one gauge field is randomly selected and updated according to the transition probability. The gauge field is initialized with state |0 everywhere and warmed up without measurements for a fixed number of iterations. After the warm-up phase, each iteration includes a measurement of the observables.
The electric energy is not diagonal in the gauge field basis. Instead of evaluating the full electric energy, we focus on the expectation value P . P acts as a lowering operator on the gauge field states. Thus, we have to evaluate an expression that has a modified gauge field on one link. We can transfer that modification to the covariance matrices by evaluating the integrals in Grassmann variables directly. The estimator for P in a Z 3 gauge theory is where Γ in is a modified version of Γ in that differs from the original one on link . Details about the calculation are provided in Appendix D.

B. Evaluation of gradients
The evaluation of gradients with respect to the parameters in T enables the efficient minimization of observables. Instead of directly tracking the derivative of the parameters through the state construction, we derive the matrix equations obtained for the covariance matrices with respect to the variational parameters. The covariance matrix of the fiducial state D does not change during the Monte Carlo computation and is the only one that contains variational parameters α ∈ {y, z}. Thus, we can calculate the gradient for an arbitrary observable O whose estimator F O (D) may depend on the covariance matrix D of the fiducial operator explicitly Since we are interested in finding the best ground state approximation with our ansatz, we calculate the gradients of the energy. They consist of two parts, the gradient of the magnetic and the gradient of the electric energy. In the case of the magnetic energy, the first term on the right-hand side of (27) vanishes since the gauge field has no explicit dependence on the parameters. It remains to calculate the expression ∂ ∂α |Ψ(G)| 2 since we know the form of |Ψ(G)| 2 from the evaluation of the transition probability (24) already. Using Jacobi's formula we obtain Combining (25) and (29), we find where ∂D ∂α is the explicit derivative of the covariance matrix of the virtual modes with respect to parameter α. This expression can be derived analytically.
In contrast to the magnetic energy, the electric energy depends explicitly on the parameters of the ansatz. Thus, the first term on the right-hand side of (27) does not vanish. The explicit form of the gradient is stated in Appendix D.
where ξ(i) is the weight for the gradient in dependence of the step. We used ξ(i) = 0.01 · 0.99 i in our simulation. The choice of parameters and the schedule of ξ(i) may be further optimised.

V. RESULTS
Applying the ansatz developed in Ref. [24] to a physical Hamiltonian, we want to ensure that we are able to capture relevant physics despite the small number of parameters of the states. In particular, we want to demonstrate that a higher number of layers leads to an improved expressibility.
As a first step, we compare to a small system with L = 2, i.e. 4 plaquettes, that can be solved with exact diagonalization (cmp. Figure 5). Due to the small lattice size, we can contract the GGPEPS exactly and do not have to use Monte Carlo. The figure and the inset show good agreement for states at high couplings where the electric energy is the dominant contribution in the Hamiltonian (1). The ground state of the electric Hamiltonian is the state with no electric excitations, i.e. the electric field is zero on all links. We expect to approximate it well because it is the state that we obtain if the operator A is equal to the identity. This happens if both parameters y = z = 0: T (y = 0, z = 0) = 1. We ob- served that the values of y and z approach zero as the coupling increases. While the high coupling regime matches well to the exact values, the low coupling regime, which is dominated by the magnetic energy, is more challenging. States with few layers show a divergent behavior at low couplings. The quadratic divergence is caused by a lack of expressibility of states with few layers: The parameters approach a constant for low coupling and the 1/g 2 term in the Hamiltonian leads to the divergence. An increase in the number of layers helps to systematically improve the states while only linearly affecting the runtime.
The error around the transition g ≈ 1 does not decrease when additional layers are used. We attribute this behavior to the specific ansatz that we are using. We do not expect a Gaussian PEPS based ansatz to hold at criticality. Figure 6 shows the energy density of the system for different lattice sizes for three layers of the parameters. Due to the larger system sizes, we cannot contract the GGPEPS exactly. The Monte Carlo estimation uses 10 4 steps for the warmup phase that is performed without measurement and 10 5 steps for the sampling. Since the Monte Carlo has to be performed for each variational minimization step, the number of Monte Carlo steps with measurements is kept rather small. Especially the calculation of the electric energy, which features a Pfaffian, is expensive.
The estimates agree very well with the ED data for an L = 2 system over a large range of the coupling. The deviations at the phase transition due to the ansatz as described above. The deviation at very low coupling for large system sizes originate from the fact that the minimization becomes increasingly costly. Especially the calculation of the Pfaffian in the electric energy is computationally expensive. While all determinants that appear in the calculation of norms can be calculated by updating previous results if the gauge field is changed, the Pfaffian has to be recalculated in every step. The Pfaffian is the single most expensive step in the algorithm. Since we are plotting the energy density in relation to a L = 2 system, deviations can be either finite size effects (in which case the MC points would be more correct than ED) or errors due to the Monte Carlo sampling procedure.
Following previous works, we expect the theory to have two phases [28,34]. According to Elitzur's theorem [39], the expectation value of any operator that is not gauge invariant will vanish, and thus a local order parameter is ruled out. Instead, following Wegner and Wilson [9,40], we can analyze the correlation in the different phases by studying the Wilson loop. The corresponding operator is gauge invariant and shows different scaling in the different phases of Z N theories. In the low-coupling regime, which is dominated by the magnetic part H B of the Hamiltonian, the expectation value of the Wilson loop follows a perimeter law which, to lowest order in perturbation theory [34], reads Here, κ p is a constant and 2(R 1 + R 2 ) is the perimeter of the Wilson loop. The scaling changes in the high coupling regime, where the electric energy is the dominant contribution to the total energy and the Wilson loop operator scales with the area of the curve. The area scaling reads to lowest order in perturbation theory [34] W where σ is the string tension. Since the potential of static charges, i.e. charges that are not dynamically coupled to the gauge fields in the Hamiltonian, increases linearly with the distance in this phase, it costs an infinite amount of energy to separate two static charges. The two static charges are confined. We can use the states that we obtained using the VMC procedure for an L = 6 lattice to evaluate the scaling behavior in the different regimes (cmp. Figure 7). As before, we used three layers in the minimization. The Wilson loop expectation values are re-computed for the minimal parameters with 10 4 warmup steps and 10 6 sampling steps. By fitting (33) to different Wilson loops W (R 1 , R 2 ) of a maximal size of L/2 and |R 1 − R 2 | < 1, we can obtain the string tension of the states. The result of the fits for different couplings is shown in Figure 7. The vertical orange line displays the expectation for the phase transition. The Z 3 gauge theory can be mapped to a three state Potts model [29] and the first order phase transition has been studied with Monte Carlo [30]. The plot shows that the string tension is almost zero in the low-coupling phase and rises to a finite value in the highcoupling, confining phase. Around the transition region, the minimization becomes difficult due to the Ansatz we are using. Thus, results in direct vicinity to the transition region might not be obtained for the ground state and one has to be careful to use them for an interpretation of confining or non-confining behavior [41]. The

VI. CONCLUSION
We show that GGPEPS are promising ansatz states for Z N lattice gauge theories in two spatial dimensions. Since the transition probability between two configurations of the gauge field is given by the squared norm of a state, the sign problem is avoided. The norm as well as the gradients for a given set of parameters can be efficiently computed with the covariance matrix formalism leading to a scalable algorithm.
By contracting small systems exactly we show that the states themselves capture the relevant physics well although they are based only on a small number of parameters. We demonstrate a systematic improvement of the energy by increasing the number of virtual fermions on the links while impacting the runtime only linearly.
The variational optimization with Monte Carlo is very successful for large couplings, but gets increasingly difficult for smaller couplings and larger lattices. In this regime, the states have to approximate states dominated by the magnetic interaction in the Hamiltonian. Since the ansatz is based on the electric vacuum on the links, this regime is challenging. Additionally, larger lattices lead to higher runtimes, especially in the calculation of the Pfaffian in the electric energy.
We expect to be able to improve the results of the Monte Carlo simulation further by changing to a more advanced sampling scheme. Currently, the algorithm updates only one spin at a time, which leads to a smaller relative change if the system size increases. The usage of collective cluster updates [42,43] or hybrid Monte Carlo techniques [11] may lead to better convergence.
Additionally, the ansatz introduced in Ref. [24] allows for static charges and dynamic fermions. The introduction of static charges allows to measure the string tension directly as an observable between two opposite charges and leads to another measure of confinement which is especially beneficial at large couplings. Simulating dynamic fermions presents the interesting possibility to study the behavior of mesonic strings.
Finally, the optimization in the weak coupling regime could be improved by starting from a different initial state on the links. If the state on the links is more suited for the magnetic Hamiltonian, the physics of the magnetic phase might be easier to capture with fewer layers.
Appendix A: Derivation of T The fiducial operator (12) used in the GGPEPS construction (9) determines the symmetries of the state |Ψ . We demand rotational invariance by π/2, translational invariance when shifting by two sites due to the staggering and charge conjugation invariance if we shift by one site. Since the parametrization was originally developed to accommodate a U (1) gauge theory [21], the formulation obeys, additionally, a global U (1) symmetry. Here, we state only the result with y, z ∈ C. y and z are the only two independent parameters that remain. The matrix is given in the mode order {l, r, u, d}. The rows correspond to the modes {l + , r − , u − , d + }, and the columns to {l − , r + , u + , d − }. In this work, we restrict ourselves to y, z ∈ R.
Appendix B: Formalism with multiple layers We achieve a higher expressibility of the ansatz states by increasing the number of virtual fermions on the links. Different layers of virtual fermions do not interact with each other and have independent sets of parameters y (i) and z (i) , where i is the index of the layer. They can be seen as different PEPS coupled to the same gauge field. Thus, the norm of the state |Ψ is the product of the norms of its layers |Ψ i where i is the index of the layer and runs from 1 to the number of layers. This construction leads to a linear scaling with the bond dimension. The matrix size of the covariance matrices stays unchanged because we do not add the parameters to the T matrix. Instead, we consider multiple covariance matrices generated by different matrices T i . Thus, we have to perform parts of the calculation multiple times with varying covariance matrices of the same size.
Since we layer only the virtual fermions, the computation of diagonal observables in the gauge field does not change. Observables like the electric energy, however, need more consideration. Due to the product structure of the ansatz state, we can write the estimator of the electric energy as a product F el = i F Finally, the gradients for the squared norm and the explicit derivative of the electric energy have to be adapted. The derivative of the squared norm enters the equations only as a fraction of the squared norm (cmp. equation (27)), we only have to adapt the expression Here, we move the derivative with respect to parameter α i ∈ {y, z} of layer i to the respective layer i since all other parameters are independent of α i . The gradient of the electric energy is adapted in a similar fashion because the derivative acts only on one of the layers.
to calculate the squared norm of the state with equation (C3). During one Monte Carlo run, D stays constant and can be calculated during the initialization. Changing the gauge field value on a link only alters Γ in (G). We refer to Ref. [24] for more details on the Gaussian mapping.
In order to calculate the squared norm of the wave function, we use the following identities [36]: where M is a complex antisymmetric 2n × 2n matrix and [X] G,θ is the Grassmann representation of the operator X in terms of Grassmann variables θ. Equation (25) follows directly from (C3).
Appendix D: Calculation of the electric energy and its gradient for ZN

Calculation of the expectation value of the electric energy
Since the electric energy is not diagonal in group element basis, we cannot use the equivalent of (23) directly. Due to the translational invariance of the states and the Hamiltonian, it is sufficient to calculate the expectation value of the electric energy over one link . The notation for Ψ(G) introduced in (9) is changed to distinguish between the group element q on link and all other group elements G to Ψ(q, G). In the following, we focus on the calculation of the expectation value P ; the extension to P + P † which appears in the Hamiltonian(1) follows directly. Since we are only considering a single, fixed link for the rest of the calculation, we drop the index .
where F el (G) is the Monte Carlo estimator of the electric energy. From the second line to the third line we use that P acts as a lowering operator on the gauge field states. The remaining expression is the product of two wave functions that differ in terms of the gauge field on one link. Using the explicit formulation of the state, we obtain (product symbols as in (9)) Ψ * (G, q )Ψ(G, q) Thus, we calculate the expectation value of the new operator U (q) ω with the density matrix resulting from the original wave function Ψ(G). Since we gauge only the right and upper modes, we can focus on the gauging transformation U (q ) = exp iΦr † + r + = exp iΦr † r with Φ = ±δ. Without loss of generality, we choose a right mode for the computation. We consider only positive modes r + for simplicity. The negative modes r − are gauged with the same expression where Φ is substituted by −Φ. For increased readability, we will skip the plus and minus signs of the modes in the following calculation.
where t = tan Φ 2 . The covariance matrix M (Φ) in (D5) of the r and l modes replaces a part of the original covariance matrix Γ in that belongs to the link that U (q) acts on. Since one link consists of positive and negative modes, we will have to substitute the single link with the direct sum M (Φ) ⊕ M (−Φ).
Due to the modification of the original covariance matrix for the projectors, we have to adapt the calculation for the overlap of two wave functions. While the identities (C3) still hold, formula (25) cannot be used. Instead we calculate the overlap using which follows from (C3). Here, X and Y are operators and Γ X and Γ Y are the covariance matrices of X and Y in terms of Grassmann variables. If the operators are Gaussian, these representations coincide with the covariance matrices in terms of Majorana fermions.
The Grassmann representation of the involved operators are Here, Γ in ( ) is the covariance matrix of link . Thus, we have to use an adapted prefactor for (D6): Tr U † q ωρ = 1 2 (1 + cos(Φ))2 −n Pf (D) Pf Γ in − D −1 , where Γ in is the modified covariance matrix of the links as defined in (D8). In the case of a Z 3 gauge, we know that cos(Φ) = − 1 2 and obtain This expression can be further simplified since the Monte Carlo estimator (D1) divides by the square of the norm and we obtain This is the expression stated in the main text as equation (26). In the case of a pure gauge theory, (D10) can be further simplified with D −1 = −D.

Calculation of the gradient of the electric energy
In contrast to the calculation of the gradient of the Wilson loop, we cannot neglect the first term in (27). The estimator of the electric energy depends explicitly on the parameters of the ansatz. Thus, we have to build the derivative of F el (D10), the estimator of the electric energy, with respect to the parameters α ∈ {y, z}.
As above, the expression for ∂D ∂α is an analytical expression. Since D is a covariance matrix of Majorana fermions in a pure gauge theory, D −1 = D † = −D holds. Thus, the first trace of (D11) is zero.