Finding the ground state of a lattice gauge theory with fermionic tensor networks: a $2+1d$ $\mathbb{Z}_2$ demonstration

Tensor network states, and in particular Projected Entangled Pair States (PEPS) have been a strong ansatz for the variational study of complicated quantum many-body systems, thanks to their built-in entanglement entropy area law. In this work, we use a special kind of PEPS - Gauged Gaussian Fermionic PEPS (GGFPEPS) to find the ground state of $2+1d$ dimensional pure $\mathbb{Z}_2$ lattice gauge theories for a wide range of coupling constants. We do so by combining PEPS methods with Monte-Carlo computations, allowing for efficient contraction of the PEPS and computation of correlation functions. Previously, such numerical computations involved the calculation of the Pfaffian of a matrix scaling with the system size, forming a severe bottleneck; in this work we show how to overcome this problem. This paves the way for applying the method we propose and benchmark here to other gauge groups, higher dimensions, and models with fermionic matter, in an efficient, sign-problem-free way.


I. INTRODUCTION
Tensor network states (TNS), and matrix product states (MPS) in particular, have been very fruitful in dealing with complicated quantum many body models in condensed matter physics. MPS are ansatz states with a built-in area law for the entanglement entropy [1][2][3] equipped with algorithms which scale polynomially with the system size, rather than exponentially. This opens the way for a very efficient ground state search [4,5], as well as studying the dynamics [6] and thermal states [7,8] of quantum many-body systems. The idea of tensor networks can be extended to higher spatial dimensions, when MPS are generalized to projected entangled pair states (PEPS) [2,3,9,10]. However, in spite of the great numerical success of MPS, PEPS in higher dimensions are generally very difficult to handle (i.e. contract) numerically; while they have been successfully applied to two dimensional models in some cases (see, e.g. [11,12]), the computational time required for contracting PEPS generally scales unfavorably when increasing the spatial dimension to more than one [13].
Gauge theories are at the heart of our modern physics. They play a central role in the standard model of particle physics, giving a local description of the fundamental interactions [14]. In condensed matter physics, they offer effective and emergent descriptions of intriguing many-body phenomena [15]. The local symmetries upon which they are based allow for a simple and local formulation of complicated interactions, and give rise to a highly constrained Hilbert space. Gauge theories exhibit some fascinating features, such as running coupling, which can be both useful and challenging: for example, in Quantum Chromodynamics (QCD), the theory of the strong nuclear force, the high energy limit is asymptotically free [16], allowing one to study collider physics using perturbation theory. The low-energy side, on the other hand, is strongly interacting and highly nonperturbative, giving rise to beautiful yet challenging physical phenomena such as quark confinement [17,18] involving many open puzzles.
A very successful method to address nonperturbative gauge theory physics has been within the framework of lattice gauge theories [17,[19][20][21]. Based on Monte-Carlo simulations of Wick-rotated path integrals in Euclidean spacetime [22,23], many challenging computations of static quantities in gauge theories have been carried out (for example, much of our knowledge of the hadronic spectrum these days is the result of lattice Monte-Carlo computations [24]). On the other hand, these methods do not allow us to directly study real-time evolution, and thus are not suitable for dynamical analysis. Additionally, in many important physical scenarios with a finite density of fermionic matter, the sign problem [25] makes the use of Monte-Carlo methods impossible.
In the last decade, some quantum information based methods have been proposed to deal with these issues. First, quantum simulation [26] offers one approach to simulate this challenging area of physics using atomic, molecular, optical and solid state devices; one can design [27][28][29][30][31][32][33] and build [34][35][36][37][38][39][40][41] quantum simulators of lattice gauge theories which are free of the above-mentioned issues, and this is a very promising and exciting avenue of research. Tensor network states, on which this work focuses, provide another quantum-information-based way to deal with them [42,43].

arXiv:2211.00023v2 [quant-ph] 11 Jan 2023
In a single space dimension, Matrix product states, combined with DMRG (density matrix renormalization group) techniques [4,44] have been successfully applied to lattice gauge theories, both in particle and condensed matter physics (see reviews [29,30,42,45] and references therein). A reliable but limited approach is to extend MPS-based algorithms to two-dimensional systems defined on ladders or cylinders. While this approach has been successfully applied to study gapped phases of matter and the transitions between them (e.g., [46][47][48][49] and others), the computational cost of this approach grows exponentially with additional dimensions, making the extrapolation to the thermodynamic limit generally impossible. In higher dimensions, tensor network techniques have been used for studying lattice gauge theories. Lattice gauge theories have been studied numerically using various methods in [43,[50][51][52][53][54][55].
In parallel, some analytical work has been done out on the properties of gauge invariant PEPS, and in particular on gauging mechanisms [56,57] which allow one to minimally couple gauge fields to matter-only PEPS in a straight-forward manner even for arbitrary groups and space dimensions. The latter method, when combined with Gaussian fermionic PEPS [58], was used for introducing Gauged Gaussian Fermionic PEPS (GGF-PEPS) [59,60]. In such states, a free (Gaussian) fermionic matter state with a global symmetry is gauged in a way analogous to conventional Hamiltonian-level minimal coupling techniques [61], resulting in a physically relevant state which describes quantum matter and gauge fields interacting in a gauge-invariant way. Moreover, it was shown that this special class of PEPS allows one to contract the states and compute correlation functions efficiently using Monte-Carlo sampling [62]. For this method, the sampling probability is shown to always depend on the norm of a quantum state, and hence it is sign-problem-free. Therefore one can use GGFPEPS as ansatz states for variational Monte-Carlo [63,64] ground state search of lattice gauge theory Hamiltonians, overcoming the common sign problem experienced by lattice gauge theories, as well as the difficulty of contracting PEPS in more than a single space dimension. Previous works have variations of tensor network states as variational Ansatz for Monte Carlo, e.g. [65][66][67]. Here, we construct anstatz states that are locally gauge invariant and adapted to the gauge group in question.
In previous work [68], this method was used to find the ground state of a pure Z 3 lattice gauge theory. While successful for most values of the coupling constant, one important issue with the algorithm was the need for computation of Pfaffians [69] of matrices whose dimensions scales with the size of the physical system. This scaling posed a serious bottleneck on using GGFPEPS as variational ansatz states [70]. In this work we present a simple way to overcome this problem, and demonstrate its use for finding the ground state of a Z 2 pure gauge theory [71][72][73][74][75]. While this computational problem on its own does not suffer from the sign problem, the tech-niques demonstrated in this work can be generalized in a straight-forward way to cases which do suffer from it (e.g., coupling the very same theory to physical fermions with an odd number of flavors and imposing the gauge constraint [48,76]). This paper is structured as follows: In Secs. II and III, we introduce Z N gauge theories with a special focus on Z 2 and constructing the GGFPEPS ansatz state. Subsequently, the GGFPEPS are minimized with the algorithm described in Sec. IV. Numerical results are presented in Sec. V and we conclude in Sec. VI.

II. PHYSICAL SYSTEM
We focus on pure gauge Z N lattice gauge theories in the 2 + 1d Hamiltonian framework; i.e. the physical system is a two dimensional spatial lattice, with gauge field degrees of freedom occupying its links. Due to the absence of dynamical matter, no physical degrees of freedom are associated with the sites.
Each link = (x, i) (labelled by the starting site x and a direction i ∈ {1, 2}) hosts an N dimensional Hilbert space. On each link we introduce two operators, which satisfy the following conditions -N th roots of unity, unitarity, and Z N algebra respectively [73]. If we label the eigenstates of P by {|p } we can immediately deduce that Q is a unitary and periodic (|N = |0 ) raising operator, Similarly, one can show that the eigenstates of Q, are unitarily and periodically lowered by P , The dynamics is given by the Hamiltonian [73] H = λ 2 2 − (P + P † ) where the second sum is over plaquettes p -unit squares of the lattice -and p 1 , p 2 , p 3 , p 4 refer to the four links around p (see Figure 1 for an illustration of the notation). The Hamiltonian is gauge invariant; that it, it is invariant under local unitary transformations of the form (whereê i is a lattice vector in the i ∈ {1, 2} direction; see Figure 2). Since [H, V (x)] = [V (x) , V (y)] = 0 for any lattice sites x, y, the Hilbert space is decomposed into a set of superselection sectors labelled by eigenvalues of all V (x) operators. We shall focus here on the sector whose states |ψ satisfy The N → ∞ limit reproduces compact QED -the Kogut-Susskind Hamiltonian [19] of a U (1) lattice gauge theory [17,21]. With this analogy in mind, we refer to the first (link) term of (6) as the electric energy and to the second (plaquette) term as the magnetic one.
In this work, we focus on the N = 2 case -a pure Z 2 lattice gauge theory [71,77]. There, we can make the choices Q = Q † = σ x and P = P † = σ z , and identify . The Hamiltonian (6) then simplifies to (again following the conventions of Fig. 1). Traditionally, a link with |p = 0 is considered as one that does not carry any electric flux while a link with |p = 1 carries it. Z 2 gauge invariance, in the sector defined by Eq. 8, implies that we only discuss superposi-tions of products of P eigenstates, involving only closed flux loops.
To gain additional insight into the physics of the model, let us consider the extreme coupling limits. First, the strong coupling limit, where λ 1, for which it will be convenient to redefine the Hamiltonian asH = λ −1 H, and in which the electric term dominates whiile the magnetic term is a small perturbation. For λ → ∞, the ground state takes the form -a product state of |p = 0 on all the links. If we decrease λ a little bit, such that the condition λ 1 is still satisfied, using perturbation theory one obtains that the ground state becomes (11) The first order correction is a superposition of all the P eigenstates which contain a single plaquette's flux loop; the second order correction will be a superposition of all the possible two-plaquette excitations, and so on. Consider the Wilson loop operator [17], which in the Z 2 case takes the form where C is some close curve on the lattice, usually chosen to be rectangular. For a rectangular Wilson loop with dimensions R 1 × R 2 , the leading order of the expectation value of the Wilson loop in the strong coupling regime will be given by the A = R 1 R 2 -th order of the perturbative series of Eq. (11), and hence the Wilson loop will decay with an area law, -a manifestation of a confining phase [17,77].
On the other hand, if we take the other extreme limit of λ 1, the magnetic energy dominates. In the extreme case of λ = 0 the ground state is the magnetic one, |ψ B , also known as one of the toric code ground states [78]. It takes the form An equal superposition of all the P states satisfies the gauge invariance condition (8). Raising λ but staying in the λ 1 regime, we can again use perturbation theory to construct the state, and a straightforward calculation shows that -a perimeter law decay, manifesting a deconfined phase [17,77].

III. CONSTRUCTION OF AN ANSATZ STATE
We would like to build an ansatz state for a variational search of the ground state of the Z 2 Hamiltonian of Eq. (9), for different values of λ. What is required from such an ansatz state?
The state must be gauge invariant, and we would like to use a construction which will allow us to couple, in future work, with fermionic matter. We therefore use the construction of Gauged Gaussian Fermionic PEPS [59,60]. Furthermore, as was shown in [62], one can perform all the relevant calculations for such states rather efficiently using Monte-Carlo.
We use this ansatz state, with the Z 2 gauge group, but without physical matter (dynamical fermions). The auxiliary degrees of freedom which shall be used for the construction of the state will nevertheless be fermionic, to allow us to couple to fermions later on and enable the efficient Monte-Carlo computation. Additionally, we would like our family of ansatz states to contain the extreme λ 1 and λ 1 cases discussed above; ideally, the ansatz will be able to interpolate between them and approximate the ground states in the intermediate regime.
We construct our states as follows: for each site x, we introduce 4F auxiliary or virtual fermionic modes, associated with the edges of links intersecting through it -F modes in each direction. We denote their creation operators by , representing the right, up, left and down directions, respectively. We will unite the 4F virtual modes on a single vertex x under a general notation of the form a † α (x) 4F α=1 , and define the Gaussian operator on each site (in general, T αβ may be site dependent, but we make it uniform to enforce translation invariance). We define the virtual symmetry operator and note that We can decompose the symmetry operator to such that and similarly for u, l, d.
We further define the gauging operators U G (x, i) (for i ∈ {1, 2}) on the links, which multiply the fermionic creation (or annihilation) operators of virtual fermions associated with them with the physical gauge field operators Q on the respective links: Note that We also define, on each link, the Gaussian operators which couple the virtual fermions associated on both sides of a link and satisfy The coupling between neighboring sites is necessary because the states would remain a product state otherwise.
Finally, the fermionic Fock vacuum |Ω is invariant under all the virtual V (x), and the gauge field state |ψ E defined in Eq. (10) is invariant under all the physical gauge transformations V (x). With all the above definitions, we are ready to define our ansatz state, the fermionic gauged Gaussian PEPS as where U G = x,i U G (x, i). Let us examine this state carefully, to understand what it may describe. First, on the right, we begin with the physical state |ψ E , which is the no-flux, strong coupling limit ground state, as defined in Eq. (10). On top of it, we act with the gauge field operator O ≡ Ω| Note that this expression contains a part that acts on the gauge-fields that is not traced out. While all virtual fermions are traced out, the expression still acts as an operator on the gauge fields. To add physical fermions to the game, one can add their Fock vacuum on the right, and add them to the A (x) operator appropriately [62]. Here we only focus on the pure gauge case and shall now show how to tailor it to the requirements mentioned above. Note that both |ψ E and O are gauge invariant (this can be shown using the symmetry properties of Eqs. (18,22,24)) and therefore our ansatz is gauge invariant, that is, it satisfies Eq. (8).
The operator x A (x) creates a virtual fermionic state when acting on the vaccuum |Ω . It is a product state of the different sites x; on each x, depending on the parameters of T , we may excite some, none, or all the virtual modes associated with the four links around it. A leg with an odd number of excitations will be referred to as one which carries virtual flux, and one with an even number of excitations will not. On top of that, we act with the gauging operator U G , and then apply the resulting operator to the initial physical state |ψ E . The gauging operator multiplies each creation operator on the outgoing links (right and up) by the physical Q operator on the same link, and thus (since Q 2 = 1), links that carry virtual flux will now also carry a physical one. However, the state is still not gauge invariant, since when we take the product of all the sites we might well end up with open flux strings which violate the gauge symmetry; the projection onto x,i w † (x, i) prevents that possibility, and ensures that the fluxes are properly connected.
We would like to obtain an ansatz PEPS |ψ with the minimal F possible, which will include the extreme case ground states, |ψ E and |ψ B , defined in Eq. (10) and (14) respectively. We would also like the state to have rotational invariance, in the lattice sense (invariance under π/2 rotations). In Appendix A we show that this can be obtained for F = 2 (two virtual fermions on each link) with (where the ordering is r 1 , u 1 , l 1 , d 1 , r 2 , u 2 , l 2 , d 2 and all the parameters are, in general, complex) and (where η 4 = −1, and this parameterization was derived for the choice η = e iπ/4 -equivalent parameterizations for other possible values exist too). As we show in the appendix, setting all the parameters to zero gives rise to |ψ E ; setting them all to zero besides b 4 = − 1 16 produces |ψ B . A gauge invariant construction can be obtained as well with one virtual fermion per link, F = 1, and (as can be seen in the appendix) it captures the physics for large and small values of the coupling constant λ with good precision, though not exactly. We refer to this construction as the "minimal" ansatz. However, as can be seen numerically, it fails in the intermediate coupling regime. Therefore, in the main text we focus on the F = 2 (cf. Eq. (26)) approach, which we denote the "optimized" ansatz. For further details and a direct comparison of the two approaches, we refer to Appendix A.

IV. THE ALGORITHM
Having introduced and justified the ansatz physically, we would now like to recall why it is useful numerically and discuss the computation algorithm. We tackle this problem in two steps: first we show how to compute expectation values for a given set of parameters by combining GGFPEPS with Monte Carlo [62]. We then demonstrate how to adapt the parameters using variational Monte Carlo (VMC) methods.

A. Q eigenbasis formulation
Following [62], we begin by expressing our ansatz PEPS |ψ in the Q eigenbasis. We introduce the gauge field configuration states, which are simply product states of Q eigenstates on all the lattice's links: As eigenstates of all the Q operators, they satisfy Due to the orthonormality of the local |q state on each link, these states are also orthonormal: In this basis, the gauging operators U G (x, i) from Eq. (21) may be seen as controlled operations, transforming the fermionic operators based on the gauge field's Q eigenvalue, where As a result (and up to an irrelevant normalization factor -the PEPS is not normalized in any case), we can rewrite the ansatz state as The wavefunction is given by where The wavefunction ψ (Q) is nothing but an overlap of two fermionic Gaussian states, (we have used the fact that U Q = U † Q ). Fermionic Gaussian states are fully classified by the elements of their covariance matrices [69]; the covariance matrix consists of correlators of products of two fermionic operators (fermionic two-point, or Green's, functions). They are central in our algorithm, based on conventional Gaussian fermionic PEPS techniques [58].
So far, we have expressed the fermionic modes using creation and annihilation (Dirac) operators. Numerically and analytically, however, following the Gaussian formalism discussed in [69] and in the context of PEPS in [58], it is more convenient to work with Majorana modes since they yield real-valued covariance matrices. As usual, for a given Dirac mode annihilated by c and created by c † , we define two Majorana modes -that is, a system with p Dirac modes has 2p Majorana modes, labeled from 1 to 2p. The Majorana modes anticommutation relations are given by the Clifford algebra The covariance matrix of a state |Φ is given by In Appendix B we show how to compute the covariance matrix of a state expressed in terms of exponentials of creation operator bilinears acting on the Fock vacuum, as the states |ψ R and |ψ L (Q) that we work with. We denote the covariance matrix of |ψ R by D, and that of |ψ L (Q) by Γ in (Q). Since both |ψ R and |ψ L (Q) are product states, it is easy to construct their covariance matrices out of local ingredients: both covariance matrices will be block diagonal. D will be a direct sum of identical blocks, each being the covariance matrix of the state created by A (x) on a single site; Γ in (Q) is a direct sum of covariance matrices of pairs of virtual fermions on the links, which will not be identical due to the gauging.
B. The norm of the state All our computations will be based on computing expectation values of operators with respect to the ansatz state |ψ , in the form of Eq. (33). However, as this state is not normalized, we would like to show how the norm is computed. Note that -the squared norm of |ψ is nothing but a sum over overlaps between two Gaussian fermionic density matrices. This can easily computed in terms of their covariance matrices [81] giving rise to:

C. Computation of the Wilson Loop
The first expectation value we are interested in computing is that of a Wilson loop W (C) as defined in Eq. (12).
The key property here [62] is that the configuration states |Q are eigenstates of the Q operators and thus also of the Wilson loop operator: where the q ( ) values are dictated by the eigenvalues of Q operators with respect to the configuration of |Q . Therefore, where F W (C) = ∈C (−1) q( ) and we define the function Note that for any Q, 0 ≤ p (Q) ≤ 1, and that Q p (Q) = 1 and thus it is a probability density function over the gauge field configuration space. Using Metropolis sampling [82], we can compute the expectation value using Markov chain Monte Carlo sampling (MCMC) [62]. Instead of the full probability p(Q), we need only the transition probability between two gauge field configurations, The denominator of Equation (45), which is hard to compute, is avoided. In the Monte Carlo procedure, we use a single-site update, i.e. we randomly select a single site and propose a new gauge field for it. Since the changes in the covariance matrices are only local, we can use the matrix-determinant lemma and the Woodbury identity to update inverses and determinants locally after each step.
In general, the exact contraction of a general PEPS is exponentially hard [83]. Here, since we picked the subclass of gauged Gaussian fermionic PEPS, we can perform the contraction required for Wilson loop computation efficiently, using covariance matrices and Eq. (42).

D. Computing P expectation values
The next operator whose expectation value we would like to compute is P on a given link. It does not act diagonally on the configuration states |Q , and thus this has be done with caution [62].
First, note that for a given link , P changes q on the link and does not affect any other links; thus, the configuration Q is identical to Q everywhere but on , where we have the opposite q eigenvalue. If we denoteQ as the configuration of gauge fields on all the links but , such that Q = (Q, q) and Q = (Q, q − 1) (the subtraction operation is obviously modulo 2), we have (48) Therefore, if we have an efficient way to compute F P (Q), we can use Monte-Carlo techniques to evaluate the expectation value of P as well. In Appendix C we show how this can be done.
In a previous work [68], the electric energy was calculated by explicitly transforming the expression to Grassmann variables. The resulting equation contains a Pfaffian that depends on the system-size. In contrast to the determinants and inverses that are used in the algorithm, the value of the Pfaffian cannot be tracked across Monte Carlo updates. Thus, the computation of a system-sized Pfaffian is necessary with each measurement.
In this paper, we introduce a new way to compute the electric energy that only depends on Pfaffians of constant size if F > 1. For F = 1, the computation does not depend on Pfaffians at all (cf. Appendix C). Here, we use the properties of the Gaussian mapping for covariance matrices to obtain the numerical value of the electric energy. In the case of the optimized ansatz (F = 2), Pfaffians of constant size enter the computation, but since they do not scale with system size, they do not hamper the computation.

E. Looking for the ground state
We wish to find the ground state by minimizing the expectation value of the Hamiltonian given in Eq. (9). Thanks to the translational and rotational invariance of the Hamiltonian and our ansatz, we can express the energy to be minimized as where n l is the number of links in the system and n p is the number of plaquettes, and P , Q 1 Q 2 Q 3 Q 4 refer to the expectation value on one particular link and plaquette which we can choose arbitrarily. The Monte Carlo procedure described above enables us to compute these expectation values for a given set of parameters α. In addition to the evaluation, we need a minimization step that drives the parameters towards the groundstate.
By computing the gradient of the energy with respect to the parameters dE dα , we can use the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [84] to minimize the parameter values.

V. RESULTS
In section III, we introduced a gauge invariant ansatz state depending on complex parameters. Its expressive power depends on the number of virtual fermions F on the links. In the following section, we will numerically benchmark the state and explicilty compare the minimal construction (F = 1) with the optimized construction (F = 2). Following analytic arguments (cf. Appendix A), we expect the optimized ansatz to match better with exact results.
The key part of the numerical computation of the energy is the evaluation of the sum in Eq. (33). In general, the number of terms in the sum scales exponentially with the lattice size. Thus, for large systems, we cannot expect to evaluate the sum exactly and we resort to Monte Carlo (MC) computations. For small systems, however, an exact contraction (EC) of the states is feasible. The exact evaluation on small systems decouples the error that we introduce by sampling with Monte Carlo from problems with the ansatz itself: even if we evaluate a bad ansatz state perfectly with MC, it stays a bad ansatz state.
As a first step, we compare the result of an exact contraction GGFPEPS with ED data on a 2 × 2 system (cf.   3). The notation 2 × 2 corresponds to a single plaquette that is closed with periodic boundary conditions, leading to 4 plaquettes in total. Here, exact diagonalization refers to solving the time-independent Schrödinger equation explicitly by diagonalizing the Hamiltonian. This is possible for a Z 2 gauge theory since the link Hilbert spaces have finite dimension.
For high couplings, where the electric term dominates, the variationally minimized data agrees well with the ED data. This behavior is expected since the electric ground state (no magnetic term) can be exactly represented with the minimal approach. We expect that the minimal approach has larger problems in the low-coupling region (dominated by the magnetic energy). In contrast to earlier works [62,68], however, we see a good agreement at lower couplings as well. The main difference stems from allowing the parameters in T to be complex. For a more detailed analysis, we refer to Appendix A.
The lower panel of Figure 3 shows the relative error where the subscripts of the expectation values indicate the method of computation. The plot illustrates that the convergence of the optimized ansatz is not only better in the transition region, but across the entire coupling region. Since we know from the analytic considerations above that the groundstate of the magnetic term cannot be exactly represented by the minimal approach, we take a more detailed look at the different terms of the Hamiltonian. The minimization of the total energy of a system is easier than obtaining the correct results for other observables. Figure 4 shows exact contraction data of the two approaches for a system of size 2 × 2. To further study the problem of the minimal approach, we visualize the total energy H, the electric H E and the magnetic H B with different colors. The curves of the total energy are the same as in the upper panel of Fig. 3. The minimal approach matches decently for the low coupling region and well for the high coupling region. In the transition region around g = 0.8, however, the total energy is far from optimal and the decomposition into electric energy and magnetic energy does not follow the actual groundstate (given by ED) at all. Note that the gray points (electric energy) and the olive points (magnetic energy) are allowed to lie under the exact solution (solid line). The variational principle holds only for the total energy and not for individual parts of the total energy.
The data for the optimized approach shows two distinct advantages over the minimal approach. Firstly, it fits much closer to the exact solution, especially in the transition regions. Secondly, the variational results follow the magnetic energy and the electric energy much better. Thus, the optimized approach is a more faithful description of the actual ground state.
The exponential scaling of the number of configurations renders the exact contraction scheme for larger systems extremely computationally expensive. For these systems, we use Monte Carlo sampling.
All variational Monte Carlo data shown in the plots is obtained with 10 5 warm-up steps and 10 5 measurement steps. Our ansatz is translationally invariant and the number of parameters scales only with F and not with the system size L 2 where L is the linear extent of the lattice. Thus, we can use the results of EC computations as starting points for the Monte Carlo computations of larger system sizes. The error bars on the variational Monte Carlo data are computed with a re-binning analysis to take into account the finite auto-correlation from the single-site update.
In the previous paragraphs, we established that the ansatz is well-suited to describe the physics of Z 2 lattice gauge theories. In Fig. 5, we check the variational Monte Carlo sampling procedure. We compare data for a 4 × 4 system computed with Monte Carlo sampling for the GGFPEPS and with an exact diagonalization code exploiting the symmetries of the system [48], which relies on the ED library QuSpin [85,86]. The data shows good agreement of the Monte Carlo sampling procedure with ED data. Finally, we can start extending the system to sizes that are not accessible with exact contractions or exact diagonalization. In Fig. 6, we show the energy of the optimized ansatz (F = 2) for different sizes. For better comparison, the energy is scaled to the number of vertices L 2 on the lattice. Additionally, we show data of standard quantum Monte Carlo (QMC) on a 6 × 6 as a comparison [76]. The good agreement between the L = 4 and L = 6 data can indicate that finite size effects do not have a large effect on the energy of the system. The situation might be different for other observables.

VI. SUMMARY AND CONCLUSIONS
In this paper, we constructed an efficient variational ansatz for the pure Z 2 lattice gauge Hamiltonian. The ansatz in form of a GGFPEPS explicitly fulfills the gauge invariance of the system. The expressibility of the ansatz can be controlled by the number of virtual fermions on the links. We show analytically that a single virtual fermion on the links (F = 1) is not sufficient to represent the weak-coupling groundstate of the theory.
Numerically, we check the analytical result by exactly contracting the state and showing that the optimized approach represents the energy more faithfully. Due to the Gaussian character of the state, the contraction is efficient and larger systems can be handled via Monte Carlo sampling.
Additionally, we demonstrate a new method to compute the electric energy for GGFPEPS. In previous work [68], the computation demanded the computation of a system-sized Pfaffian in every measurement. This computation is substituted by a matrix multiplication and an inversion which can be tracked through the local sampling procedure. This can enable the exploration of larger systems and removes one of the big runtime penalties of the algorithm.
One of the immediate next steps to add physical fermions to the system. The formulation of the state is written in terms of fermions to allow a seamless integration of physical fermions in the ansatz.
Furthermore, the more efficient formulation algorithm introduce here, may enable simulations in three space dimensions, as well as of other gauge groups, including compact ones such as U (1), SU (2) or SU (3), and quite possibly serve as a way to study non-perturbative gauge theories, in particular 3 + 1d QCD with a finite chemical potential, which suffers from the sign problem.

Rotational Invariance
First, we would like to ensure that the state is rotationally invariant, in the lattice sense, that is, invariant under π/2 rotations. We will do so by briefly reviewing the procedures of refs [59,60]. Let us first define how rotations are carried out. Given a lattice site, x = (x 1 , x 2 ), we define its rotation by We can then define the rotation of physical operators simply as where R p is the unitary operator implementing the π/2 rotation of physical operators and states, and similar relations hold for the P operators. Clearly, the strong coupling vacuum is rotation invariant, and therefore it is easy to see that the weak coupling vacuum |ψ B is invariant too, following its definition in Eq. (14). We also define the rotation of the virtual degrees of freedom, implemented by the unitary R v , by where R αβ is a matrix which relates to the rotation of the legs: taking r to u, u to l etc. Clearly, it has to be unitary. Thus, it has to be a permutation matrix (up to phases). If we choose the ordering and assume that the rotation does not create any mode mixing, we can simply write R as the direct sum where where |η| = 1. If we wish to couple to physical fermions, the virtual fermions should have similar transformation rules [59]; since a complete rotation puts a minus sign on a fermion, this phase must satisfy Note that the virtual vacuum is invariant under rotations, and thus we can guarantee, following Refs. [59,60], that our PEPS |ψ , as defined in Eq. (25) will be rotationally invariant if as well as The rotation property of A (x) from Eq. (A9) is obtained if and only if the matrix T satisfies the equation (just by acting with the rotation operators on the exponential of A (x) explicitly). This further constrains the parameterization. The rotation rules of w (x, i) may be demanded in a similar fashion. Using the definitions from Eq. 23, we get that the proper transformation rules (A10) are obtained if and only if giving rise the consistency equation Since we do not wish η to depend on F , let us consider what happens for F = 1 where W 1,2 are simply numbers and their transposition is irrelevant; this forces us to satisfy Eq. (A7) even in the absence of physical fermions.
We choose not to include any free parameters in W 1,2 ; a common trick in PEPS theory [3] allows one to absorb all the free parameters into the on-site tensors (T in our case).
With all that at hand, we can now proceed to construct the most suitable ansatz state |ψ while meeting the required contstraints.

Minimal ansatz
The most minimal construction we may try is F = 1: a single virtual fermion per leg. Then, T αβ is a four dimensional complex matrix, and a † α has four components, one on each leg (r, u, l, d). Due to the fermionic anti-commutation relations, we obtain that T is an antisymmetric matrix, reducing the number of allowed complex parameters to six. Solving the rotational invariance conditions of Eqs. (A11) and (A12) reduces the number of free complex parameters in T to two: and gives rise to It is easy to see that the strong coupling vacuum is included here, simply by choosing y = z = 0, but which other states can be created with this choice? For that, we attempt to understand the meaning of the y, z parameters. Recalling that the columns and rows of T are ordered as {r, u, l, d}, we can see that the parameter y relates to the creation of virtual fermions along a straight line (r and l, or u and d) while z is associated with corners. To see this clearly, we can rewrite our A (x) operator as (we omitted the coordinate x which is shared by all operators for simplicity, and the graphical notation should be taken with caution, as it does not account for the ordering of creation operators -one can see the first row as defining proper ordering for the diagrams below). We can open the brackets, and discover, due to the fermionic statistics and the fact the fermionic modes cannot be excited twice, that we have the following possibilities: The PEPS |ψ is created by a product of such operators on all sites, acting on the Fock vacuum and the strong coupling ground state |ψ E . The gauging operator adds Q on all the links where r or u are excited, and thus makes the "physical flux map" look like the virtual one. The remaining action of w † operators and projecting back on the virtual vaccum, completing the definition of the PEPS in Eq. (25), guarantees gauge invariance and proper closing of the physical flux loops. We can use this to consider the structure of |ψ as a series in the parameters y, z (without assuming that any of them is small).
The zeroth order will be obtained by the 1 ingredient of all A and w operators, and thus it will simply be The next order will have to involve the shortest flux loop -a plaquette -requiring a contribution of four sites which is not 1, and 1 everywhere else. This will be a combination of four corners, and thus will be accompanied by 16z 4 . The final phase will be determined by the fermionic anticommutation rules when contracting with the right terms from w † ; a straightfoward computation shows that We define the second order as that with two plaquettes excited. Here, there are several possibilities; first, of plaquettes p, p with no link or site in common. Then it is straightforward to show that the amplitude is the square of that for a single plaquette: −16z 4 2 = 256z 8 . Another option is that of two plaquettes which share a link; then, we need four corners and two straight lines, and the amplitude per such a single excitement is −64z 4 y 2 ; finally, we consider plaquettes with one site in common, which can be created in two ways, giving rise to two contributions to the amplitude: eight corners (contrbuting 256z 8 ) or six corners and two straight lines (−256iz 6 y 2 ). Continuing in this way is possible but tedious, so we stop here in order to evaluate what this state is able to describe. First, assume that y = 0 and |z| 1. Then the PEPS may be seen as a perturbative expansion, and we have If we pick z = −128λ 2 −1/4 1, we get, in leading order, the perturbative solution of the ground state for λ 1 as in Eq. (11). The next order, however, cannot be obtained, since we cannot excite two neighboring plaquettes with y = 0, so we will not extend the discussion on this perturbative limit further.
Next, we consider the weak limit. We wish our ansatz to cover both extreme limits, so what about |ψ B ? We can take a look at its definition in Eq. (14) and start expanding all the brackets. We can clearly see that the superposition includes |ψ E as we have in our PEPS |ψ and that all the orders (in particular the order of a single plaquette), no matter how many plaquettes are excited and where they are, have the same amplitude 1. Going back to what we have just derived, this implies that we need z 4 = − 1 16 . Going up to the second order, there are several kinds of terms; the pairs of far plaquttes carry an amplitude of 256z 8 = 1 as required; but those who share a link carry an amplitude of −64z 4 y 2 = −4y 2 , thus we require y 2 = − 1 4 . However, we will run into contradiction with the other kind of terms, pairs of plaquettes with only a single site shared. Their amplitude will be 256 z 8 − iz 6 y 2 = 1 ± 1 = 1, and our conclusion is that the minimal ansatz does not cover the weak limit and thus does not satisfy our basic requirements from an ansatz.

Optimized ansatz
Next, we try to build an ansatz that will include |ψ B , the weak limit ground state too, as defined in (14) with two virtual fermions per leg (F = 2). This time, we shall use a rather more constructive, bottom-up approach, by first considering a single plaquette.
Consider one plaquette and label the sites around it by a, b, c, d, as in Fig. 7. On each leg of the links we introduce a single virtual fermionic mode, eight altogether, and we name them as in the figure. Note that the sites a, b, c, d are entirely independent of the parameters a, b, c, d in Eq. (26).
Suppose this is our entire system, and we construct our PEPS for it; the initial gauge field state will simply be |0 = |p = 0 ⊗4 . The A operators will take the most general forms Gauging will be done as usual. In order to connect the modes properly, and following the previous construction, we define the operators Finally, the PEPS for a single plaquette will take the form It is easy to verify that the contraction of the virtual fermions will give rise to Why did we bother to do all that? Because we wish to generate an operator which is a product of such operators on all the plaquettes, Clearly, when f p1 f p2 f p3 f p4 = 1 for all the plaquettes, In order to build this O, we consider a PEPS |ψ with F = 2 -two virtual fermions per mode. On each site x we set which allows a site to connect to all the four plaquettes around it -to one plaquette it plays the role of a from the above construction, and to the others, the role of b, c, d.
(A27) which accounts for the link being ab for the plaquette on top of it, and dc for the one beneath it. Similarly, (A28) Plugging these into the PEPS construction |ψ gives us All we have to do is to show that this is embedded in some general parametrization with F = 2. Such a construction will have an eight dimensional T matrix. Demanding the rotation invariance condition (A11), we obtain the T matrix introduced in Eq. (26). For the w operators we stick to the choices of Eqs. (A27) and (A28). They satisfy the rotation invariance conditions of Eq. (A12). Then, if we set all the parameters of the T matrix (26) to zero, other than b, we can get |ψ B ; inspecting Eq. 26 carefully and comparing it with our |ψ B construction implies that Taking the product, we obtain f a f b f c f d = −16b 4 and hence setting b 4 = − 1 16 and all the other parameters of (26) to zero, gives us |ψ B . Trivially, setting all the parameters to zero gives rise to |ψ E . Therefore, we choose the F = 2 PEPS with T given in (26) and w from (A27) and (A28), or equivalently (27), as our ansatz. It satisfies all the requirements we set for ourselves: gauge invariant, rotation invariant, and including the extreme cases |ψ E and |ψ B with a minimal number of parameters.
A big difference between the ansatz used in this work and previous work [59,68] is lifting the restriction to real parameters in T . In Fig. 8, we illustrate the differences in using real/complex-valued parameters and changing the number of fermions on the links. The mention of layers in the figure refers to the idea of enlarging the number of parameters by adding additional GGFPEPS coupled to the same gauge field. In the absence of physical fermions the contraction of the GGFPEPS can be seen as a wavefunction in the group element basis (cf. Eq. (33)). If we choose independent PEPS, we can compute all observables independently with only a linear runtime penalty in the number of layers. For more details on layers, we refer to Ref. [68].
In Fig. 8, Panel (a) shows exact contraction data for the minimal ansatz with real parameters. As in a previous work [68], we see a 1/λ divergence for low couplings. This is exactly the prefactor in front of H B . Thus, the ansatz state reaches a minimum energy for H B and fails to lower it further. With an increase in the number of layers, the problem can be mitigated, but a distinct deviation exactly at the transition region remains. Panel (b) shows a similar plot for the optimized approach. Even for the minimal number of a single layer, we get good agreement with the expected curve (solid line). Upon increasing the number of layers, we see improvements in the transition region. In panels (c) and (d) the minimal and optimized approach are shown for complex parameters, respectively. We note that the divergence at low couplings can be mitigated with just a single layer also in the F = 1 case [panel (c)] for complex parameters. The transition region, however, remains hard to access for the minimal approach. Only the choice of complex parameters and the optimized ansatz of two virtual fermions on each link (F = 2)[panel (d)] leads to a good match across the full range of couplings.

Appendix B: Computation of the covariance matrix in Dirac modes
We consider a fermionic Gaussian state of the form If M ij ∈ R, U is orthogonal, but the decomposition in (B2) involves U T and not U † for both R and C.
In both cases (M ∈ R and M ∈ C), with λ k ∈ R and λ k ≤ 0. We denote M 0 in terms of a direct sum where M 0 (k) = iλ k σ y . In this canonical basis, |ψ can be written as a product of BCS states |ψ k : and b † i = U ij a † j . The covariance matrix of |ψ in this canonical basis, Γ 0 , can be written as a direct sum: where Γ 0 (k) is the 4 × 4 covariance matrix of the BCS state |ψ k . The Dirac covariance matrix for |ψ k is given as and v(k) = λ k √ 1+λ 2 k and 1 − M 2 0 (k) = 1 + λ 2 k 1. Using the structure of the direct sum for R and Q, we obtain In the original ordering of the operators b and b † we get Here, we ordered the operators such that {b 1 , · · · , b R , b † 1 , · · · , b † R }. Although all λ k ∈ R, we keep the notation M 0 for easier notation later.
We rotate back to the a basis that we started with and define In the a basis, we obtain which evaluates to The covariance matrix in terms of Majorana modes can be computed by a linear transformation from Eq. (B18) which follows directly from the definition of the Majorana modes in Eq. (37).
The direct connection between the parametrization M and the Majorana covariance matrix enables us to compute the derivatives of the covariance matrices automatically by symbolic differentiation. Since the matrices are not modified during a Monte Carlo evaluation, we can store the numerical evaluation of the derivatives for a given set of parameters.
Appendix C: Calculation of P Here we shall elaborate and give detail on the computation of F P (Q), that we need for the computation of P and hence for the electric energy part of H . Following Eq. (48), recall that (C1) If we define Ω (x, i) as the projector onto the empty state of all the virtual fermionic modes on the link (x, i), and we can express the numerator of F P (Q) as (without loss of generality -due to the rotational symmetry -we assume that we compute the electric energy for a horizontal link). The numerator is the unnormalized expectation value of the Gaussian opera- For simplicity, let us first assume that F = 1 -the minimal ansatz introduced above. If we denote byω the product of ω (x, i) on all the links but the one where the expectation value is computed, we obtain for the numerator where r, l are the annihilation operators of the virtual fermions on both sides of the link we study (belonging to two neighboring sites). In the third line of (C5), we used the fact thatω 2 ∝ω (as a non-normalized density matrix of a pure state), as well as thatω commutes with the modes in parentheses because it acts on different links. For convenience, we define the state φ =ω |φ(Q) .
The quantity we wish to compute is the expectation value of an operator acting on the l, r modes of one particular link, with respect to this state. Since |φ(Q) is Gaussian, this can be extracted from its covariance matrix. All covariance matrices so far have been formulated in terms of Majorana modes. Thus, we express the operator whose expectation value we seek in (C5) in terms of Majorana modes: In a previous work [68], we rewrote the operator in (C7) in terms of Grassman variables and obtained a new form for the gauged covariance matrix of the projectors Γ in . The resulting expression required the computation of a system-sized Pfaffian which is computationally expensive since it must be recomputed with every measurement. To our knowledge, there is no algorithm to infer the new value of a Pfaffian after a local change in a matrix [70]. In the following, we present a different method that explicitly uses the Gaussian character of the GGF-PEPS.
Using standard fermionic Gaussian states [69], and in particular fermionic Gaussian PEPS, techniques [58], we can compute the covariance matrix of |φ(Q) using the so-called Gaussian map; such a Gaussian map takes as its input the stateω, involving the modes on all the links but the one we are interested in, and it is parameterized by the covariance matrix D of the state |φ(Q) . We sort the fermionic modes into two groups: the ones which are contracted -those on all the links but the one we look at, and the non-contracted, or open ones, which are all the rest (mathematically, this is equivalent to virtual and physical modes respectively, in conventional fermionic PEPS constructions [58]) -see Fig. 9 for graphical explanation.
The output of the Gaussian map is the covariance matrix of φ , which belongs to the Hilbert space of the non-contracted fermions [87]. It is given by [58,69] where D oo , D oc andD are the blocks of the covariance D containing correlations of open modes with themselves, open with contracted, and contracted modes with themselves, respectively, as depicted in Fig. 10. Γ in (Q) is the covariance matrix ofω. Using the basis l (1) , l (2) , r (1) , r (2) for the modes, we can write the matrix Γ v as with · = φ · φ . Since Γ v is anti-symmetric, we do not fill out the lower half of the matrix. We identify the expression (C7) in terms of elements of the covariance matrix The full expression of P reads Finally, we can directly remove the anti-hermitian part since P is a Hermitian operator in Z 2 and it will add up to zero in the end. Thus Here, the norm φ|φ can be computed using (42) for the modified matricesD and Γ in .
In the case of the optimized approach with F = 2, the expression corresponding to (C7) is The subscript indices denote the different values of F in the system. While we can write the expectation value in terms of a multiplication Majorana modes, we cannot directly transform the products into elements of the covariance matrix. The products of four Majorana modes are identified with Pfaffians in terms of submatrices of the covariance matrix D using equation (17) from [69] The equation [69] tr(ρi p c a1 c a2 · · · c 2p ) = Pf(Γ v | a1,··· ,a2p ) , with 1 ≤ a 1 < · · · < a 2p ≤ 2n, enables the transformation from products of Majorana modes to matrix elements of the covariance matrix. Here, Γ v | a1,··· ,a2p is the 2p × 2p submatrix with the indicated rows and columns. More concretely, we need the contraction of 4 Majorana modes In terms of Pfaffians, we can write Eq. (C13) as The final expression of P reads, after omitting all imaginary terms, P = ψ| P |ψ ψ|ψ φ φ ψ(Q, q) ψ(Q, q) p(Q, q) .

Appendix D: Calculation of derivatives
The idea of a variational algorithm is to minimize the energy by (iteratively) adapting the parameters of a state. While gradient free optimizations are possible, the gradient of the energy with respect to the parameters, speeds up the minimization significantly. The total energy of the system consists of two parts: H = H E + H B . Since the electric energy is non-diagonal in the group element basis (cf. sec. IV D), we go through the relevant calculations in detail. The easier case of the magnetic energy (diagonal in the group element basis) follows directly.
We consider an observable O and its expectation value Our aim is to calculate the derivative with respect to a parameter α of the matrix T (α). The number of parameters depends on the value of F . The derivative of the observable can be written as The expression can be derived by considering the logarithmic derivative ∂ ∂α ln O and transforming back at the end of the calculation. If the observable does not explicitly depend on the parameters, as is the case for Wilson loops, the first term in (D2) vanishes.
In a first step, we focus on the calculation of the second and third term which are always present since the norm always depends on the parameters. Instead of directly tracking the derivative of the parameters through the state construction, we will use the chain rule. The matrix D does not change during one Monte Carlo computation with a given set of parameters and is the only one that contains the parameters.
The derivative of the norm is ∂ ∂α (D3) The calculation is described in more detail in an Appendix of Ref. [68]. According to the calculation in Appendix C, the electric energy is obtained by changing the Gaussian map. The expression for the electric energy contains the covariance matrix and two norms ψ|ψ and φ φ that depend on the parameters. The exact expression depends on the ansatz.
The derivative of the expression reads with v andṽ given by The tilde decorations used here are those introduced in Appendix C. When using the optimized ansatz (F = 2), we have a slightly different expression for the electric energy where f (Γ v ) is the sum of Pfaffians in (C16) without the prefactor of 1 16 . The derivative of the expression is structurally similar to one of the minimal approach with v andṽ defined as above.
In comparison to the minimal ansatz, the expression for the optimized ansatz contains Pfaffians of certain rows and columns of the resulting covariance matrix in ∂ ∂α f (M ). These are the results of Wick contractions from four-body correlators. Since f (M ) is a sum of independent terms, we only treat a single Pfaffian explicitly ∂Γ v | i1,i2,i3,i4 ∂α . (D9) For both, F = 1 and F = 2, the derivative of the electric energy contains the derivatives of the covariance matrix Γ v . The expression reads where we used ∂K −1 ∂α = −K −1 ∂K ∂α K −1 . Similar to the idea of tracking (D −1 − Γ in ) −1 , we can also track (D − Γ in ) −1 and avoid the expensive matrix inversions in each measurement computation.
The derivative of the matrices D oo , D oc ,D and D are derivatives of the covariance matrix of the Majorana modes D. As described in Appendix B, we know a direct construction of D in terms of T . Thus, we can calculate the derivatives symbolically before the actual computation and insert the appropriate parameters as needed.