Gauss law, minimal coupling and fermionic PEPS for lattice gauge theories

In these lecture notes, we review some recent works on Hamiltonian lattice gauge theories, that involve, in particular, tensor network methods. The results reviewed here are tailored together in a slightly different way from the one used in the contexts where they were first introduced. We look at the Gauss law from two different points of view: for the gauge field, it is a differential equation, while from the matter point of view, on the other hand, it is a simple, explicit algebraic equation. We will review and discuss what these two points of view allow and do not allow us to do, in terms of unitarily gauging a pure-matter theory and eliminating the matter from a gauge theory, and relate that to the construction of PEPS (Projected Entangled Pair States) for lattice gauge theories.


I. INTRODUCTION
Gauge symmetries, that provide the standard model's description of interactions, are fascinating. They have many interesting properties, that make them an excellent playground for endless studies of non perturbative physics. They offer puzzling phenomena to study and understand, and solve, due to their rich symmetry and non-perturbative nature; perhaps the most famous one is the confinement of quarks in Quantum Chromodynamics (QCD) and other non-Abelian gauge theories [1], but many others too.
Asymptotic freedom [2] tells us that in high energies the coupling of QCD is very small, allowing one to use perturbation theory and Feynman diagrams to describe and understand collider physics very well. However, the other, non-perturbative low energy physics side, is yet to be understood. A very successful approach to that has been lattice gauge theory [1,3], where spacetime is discretized in a way that allows one to regularize the theory in a gauge invariant manner, and to perform very successful Monte-Carlo calculations of many important physical properties -such as the hadronic spectrum [4]. On the other hand, as the Monte-Carlo calculations are done in a Euclidean spacetime, they do not allow, in general, to consider real-time evolution, and encounter the fermionic sign problem [5] in many physically relevant scenarios required, for example, for the study of the exotic phases of QCD [6].
Recently, much effort has been put into lattice gauge theories from the quantum information and optics point of view, suggesting new ways to tackle these problems. One is quantum simulation [7][8][9], which suggests to map the gauge theory degrees of freedom to those of a controllable quantum system that could be used as an analog experiment, or a digital quantum computer for gauge theories. The other involves the application of tensor network tools, such as MPS or PEPS [10,11].
Approaching lattice gauge theories in both methods is done in the Hamiltonian approach, first introduced by Kogut and Susskind [12]. Unlike the path integral methods used widely in particle physics in general, and lattice gauge theories in particular, in the Hamiltonian formulation one deals with states and operators in a Hilbert spaces, that have a very special structure due to the gauge symmetry.
Constructing a quantum simulator or a tensor network state involves a very interesting challenge: one has to reconstruct a physical model from "other", different ingredients. In this process, one faces the most fundamental elements of a theory, which is decomposed to its smallest ingredients. This, in some sense, is what we would like to discuss in this lecture: it will mostly be about enforcing gauge invariance on quantum states -the construction of tensor network states for lattice gauge theories, through a discussion of the Gauss law.
The key ingredient of gauge theories is the local nature of the symmetry: a gauge symmetry is a local symmetry, which means that they involve many (local) conservation laws, or constraints. The gauge symmetry is formulated, in both classical and quantum gauge theories, by means of the Gauss law, an equation -or a set of equations, one per space point -that relates the gauge field and the matter fields. Its simplest form, of classical electrodynamics, is -the divergence of the electric field equals the density of matter charges; it can be made more complicated when the gauge field is non-Abelian, become an eigenvalue equation for the so-called physical states after quantization, or a difference equation for lattice gauge theories, but the essence is the same in all those cases: the matter is the source of electric fields. The Gauss law is a static equation, even though its ingredients are, in general, time dependent. It involves no temporal derivatives; as the generator of a symmetry, it commutes with the Hamiltonian. It formulates constants of motion, rather than an equation of motion. Although it can be obtained from the Euler-Lagrange equations for fields, it is an artifact of the continuous Lagrangian formalism for fields, that puts space and time on equal footing. It cannot be obtained through the Hamilton equations, but rather appears as a constraint accompanied by a Lagrange multiplier once a Legendre transformation to the Hamiltonian formalism is carried out. In lattice formulations with continuous time, which do not put time and space on equal footing, it cannot be obtained as an equation of motion as well. As it is a constraint, it means we can try and solve it, and if we manage to do that, we can solve it and plug it into the other equations of motion, we can reduce the number of degrees of freedom.
In this lecture, we will discuss the nature and possibilities of solving the Gauss law, in two contexts. First, we will look at it a a differential equation for the gauge field, and see that it tells us that gauging a free matter theory (minimal coupling) could not be done, in the most general setting using a local unitary transformation. On the other hand, we will see that if we can break the system into pieces with simple Gauss laws (and we will explain what "simple" means), we can modify globally invariant Hamiltonians or states by gauging locally and unitarily each piece alone, and then tailor everything together; this will be done using the Trotter-Suzuki decomposition for Hamiltonians, and construction of PEPS -Projected Entangled Pair States. Second, we will discuss the possibility of solving it for the matter, and using that for the elimination of the matter degrees of freedom, or, in case they are fermionic -their replacement by spins.

A. Minimal Coupling and its Compact Lattice Formulation
A very conventional textbook approach to gauge theories is to consider first a non-interacting matter field theory, with some global symmetry: the reader is suggested to naively try and act on a globally invariant Lagrangian or Hamiltonian with local transformations instead of global ones, only to quickly find out that the kinetic (derivative) term is not invariant under them. Then, in order to obtain something with a local symmetry nevertheless, another degree of freedom -the gauge field -is introduced as a geometric connection, modifying the previously quadratic, non-interacting terms into interaction terms of the matter with the field, in a procedure known as minimal coupling. At this point the gauge field has no dynamics; this could be introduced by including further terms in the Hamiltonian or the Lagrangian.
The most famous example, perhaps, is the derivation of the Quantum Electrodynamics (QED) Lagrangian from that of the free Dirac theory [20]; one begins with which is symmetric under the global transformation ψ (x) → ψ (x) e iΛ . In order to make it symmetric under local transformations of the form ψ (x) → ψ (x) e iΛ(x) , one introduces the gauge field A µ (x) and replaces the regular derivative ∂ µ by the covariant one D µ = ∂ µ − iA µ (x), to obtain instead Finally the dynamics is obtained by adding the Maxwell part is the well known, gauge-invariant field strength tensor). The complete process could be summarized as One could perform a similar procedure on the lattice. Consider, as a simple example, a two dimensional spatial lattice, with continuous time, as in the usual Hamiltonian formulation of lattice gauge theories [12]. We define a fermionic mode (creation operator) ψ † (x) at each vertex (lattice site) x ∈ Z 2 . We introduce a staggered fermionic Hamiltonian, where one sublattice corresponds to particles and another to anti-particles; if one tunes the hopping coefficients in a slightly different way, which was not done here for the sake of simplicity, a (doubled) Dirac model is obtained in the continuum limit [21]: whereê i is a unit vector in the direction i = 1, 2. H f is invariant under global phase transformations, generated by the total fermionic number which is a conserved charge.
To make the symmetry local, we introduce on each link = (x, i) of the lattice a new Hilbert space, for the gauge field, with two conjugate operators: the compact, angular vector potential φ (x, i), and the electric field E (x, i) [3], E (x, i) is a U (1) angular momentum operator, conjugate to the angle φ (x, i); we recognize the Hilbert space on each link of such a compact U (1) lattice gauge theory as this of a particle moving on a ring. As phase operators are not well defined, we work instead with the group element operator -playing the role of a Wilson line along the link. On each link, the Hilbert space may be spanned, in particular, by two complete sets of states [22]: either the magnetic one, of group elements (or more accurately, in this case, parameters), |φ , with the orthogonality relation The gauge symmetry implies that -every physical state |ψ is an eigenstate of all the Gauss law operators G (x), with eigenvalues q (x) which we call static charges: This eigenvalue equation is simply the Gauss law -a quantum lattice version of the well-known continuum classical equation (1). The static charge configurations {q (x)} are constants of motion, splitting the physical Hilbert space into different superselection sectors, that are not mixed by the dynamics. When one acts with a local gauge transformation on a state, a global phase depending on the static charge configuration sector H ({q (x)}) to which it belongs will appear. Due to the superselection rule we do not discuss superpositions, and thus only global phases can appear, and it is well known that they play no role in quantum mechanics. Therefore, in some sense all the states in H are gauge invariant -but with different phases, depending on the static charges. In our context, however, a gauge invariant state will be such that is strictly invariant, i.e. without a phase. This generalizes also to other, non-Abelian groups, where static charges imply multiplets of states that are connected through gauge transformation, unless the static charges are group singlets. From now on, when we discuss the physical or the gauge invariant Hilbert space, we will refer only to the H ({q (x) = 0}) sector.
A general physical state |ψ ∈ H ({q (x) = 0}) satisfies and thus it may be expressed as a superposition of product states of gauge field and matter, with the same eigenvalues for G (x) , Q (x): However, note that once |{Q (x)} Matter is fixed, there is more than one choice for |E ({Q (x)}) Gauge (it is, in general, a superposition). Which means that a unitary gauging map of the form (minimal coupling) is, in general, not possible. On the other hand, once |E ({Q (x)}) Gauge is fixed, there is a unique choice for |{Q (x)} Matter (this is true for our U (1) staggered case, but not in general). Therefore a matter eliminating transformation, exists. This is strongly related to the Gauss law and its properties -as we shall discuss below.

C. Can Minimal Coupling be Unitary and Local?
In quantum mechanics we like unitary transformations. In particular, when two Hamiltonians are connected by a unitary transformation (just like any other pair of Hermitian operators), they will have the same spectrum, and their eigenstates are straightforwardly connected.
Let us focus on the first step, H f −→ H f . It is clear that H f , U (x, i) = 0 as no gauge field dynamics is introduced at this point. Could we then find a unitary gauging, or minimal coupling transformation U G -for which To make the situation even simpler, let us first only deal with a single link, with two fermionic modes ψ, χ on its edges. The globally invariant Hamiltonian is with a global symmetry generated by and gauging it it will lead to with the symmetry charges (Gauss law operators) Our initial Hilbert space is H f , a subspace of the fermionic Hilbert space, that contains the globally invariant fermionic states. Gauging maps it to H phys , a subspace of the product space H f × H g (where H g is the link's Hilbert space), which is invariant under the transformations generated by G ψ , G χ . We will replace the initial space by H f × {|0 }, where E |0 = 0. We would like to construct a unitary transformation U G , that maps H f × {|0 } to H phys -that is, it will entangle the vertices and link degrees of freedom, in a way that respects the symmetry generated by G ψ , G χ .
We will have to act on the state |0 in a way that will introduce the desired conservation laws. Since they have to be consistent with each other, we need some relation between the two fermionic number operator, but this is satisfied by our initial, globally invariant elements of H f . So we simply have to identify E with ψ † ψ or χ † χ − 1: that is, to raise the electric field on the link by an amount ψ † ψ, or lower it by χ † χ − 1. Let us choose the first option. This is done by the unitary [18] where |φ is a phase eigenstate, such that It is a controlled operation that could be seen either as raising the electric field E with respect to the fermionic charge ψ † ψ, or as rotating the fermionic operators ψ † , ψ with respect to the phase operator φ. While the first interpretation suits the point of view of the quantum state, the second one suits that of the Hamiltonian, explaining us intuitively that, indeed, -so we see, that unitary gauging was obtained by solving the simple Gauss law (30) which was both local and uniquely solvable.
Could we now extend this to the the Hamiltonian H f in arbitrary dimension? The answer is no. Of course, for each link = (x, i) we can define a gauging transformation such that however, each mode appears in more than one link, and applying the product U G ( ) to the whole Hamiltonian H f will give rise to a very complicated and messy expression which is not H f : since now each Gauss law involves more than one electric field, all fermionic operators will be rotated with respect to the gauge fields on all the links around them.
Let us move on more slowly then, and extend our one link system to an open line. The globally invariant Hamiltonian is with a global symmetry generated by Gauging it it will lead to (38) with the symmetry charges (Gauss law operators) Solving this equation gives us something unique, but nonlocal - so, if we like to repeat the single link procedure, we will need to put initially a trivial |0 state on each link, and then change it by an amount of y<x Q (y). This leads to a nonlocal gauging transformation, of the form that gives us what we want. Its inverse form can be used for eliminating the gauge field degrees of freedom in the one dimensional case, when one starts with the gauge invariant theory, but it will introduce some nonlocal interactions (from H E ). It can be used, for example, for the construction of quantum simulators, as was done in the trapped ions experimental realization of the Schwinger model in Innsbruck [23]. We have seen that for a single link we can gauge with a unitary -a controlled operation -because the Gauss law has a unique solution. For a one dimensional system, we can gauge by a unitary transformation again, since a unique solution exists, but we lose locality. In more dimensions we cannot do even that anymore.
In order to understand it better, and in a more general way, let us look again at the equation that is responsible for all that -the Gauss law. Either its lattice form (14), or the well known continuum form, This is the equation that manifests the local symmetry, introduced to the system in the minimal coupling procedure, relating the gauge (electric) field to the matter, that existed before. In the process of gauging, we "complete" a physical setting that only had the matter field to one that involves the gauge field too.
Suppose we wish to gauge a globally invariant Hamiltonian, or state, by introducing a gauge field that will satisfy this equation. That means that we have to solve it for the electric field. The problem is that such a solution is, in general, not unique (with the exceptions mentioned above): we have a single, non-vector differential (or difference, on the lattice) equation, for a vector field; for example, a general gauge invariant state could be a superposition of many different states satisfying the Gauss law, even if the matter charges are fixed. However, if we have no unique way to determine the gauge field configuration from the matter, it means we cannot gauge using a unitary transformation.
We see, therefore, that gauging by a local unitary will generally not work. Indeed, in the process of minimal coupling we modify the Hamiltonian, and get something which is completely new, not the result of acting with a unitary operator: we really modify the system, rather than transform it.
Hence, the only scenarios that allow one to gauge using a unitary and local transformation are those in which the Hamiltonian or the state to be gauged are constructed out of separate, non-interacting local parts, each of which could be separately gauged in a unique way. Only after gauging, one tailors the already gauged pieces together, and introduce nonlocal correlations. That is, we do not transform the whole object (Hamiltonian / state), but rather its building blocks, before connecting them. We do not transform, we modify.

D. The Case of Hamiltonians: Trotterized Gauge Invariant Dynamics
From a Hamiltonian point of view, this can be achieved using trotterized time evolution. The Trotter-Suzuki decomposition [24], which is very useful for quantum simulation [25], allows one to approximate the time evolution of a Hamiltonian by a product of Trotter steps: short time evolutions of different parts of the Hamiltonian, that do not necessarily commute. It is exact if each Trotter step is infinitesimally short; otherwise, for finitely short steps, error bounds can be computed.
For our purposes, we will break the Hamiltonians into separate link contributions, This decomposition allows us to gauge separately: such that with a finite but large N this can be used as an approximate time evolution, as widely done for quantum simulation purposes in other contexts, Why does that work? Because every Trotter step here is just like the single link Hamiltonian H 2 we considered above (27), and its gauging gives rise to simple, one dimensional Gauss laws generated by operators such as those of (30) rather than the ambiguous, vector Gauss laws of the general case (14). They are only "switched on" once a sequence of Trotter steps is considered, as symmetries of the complete time evolution.
Some errors due to the trotterization (non-commutativity of intersecting links contributions to the Hamiltonian) will arise; note, however, that these errors have nothing to do with the gauge symmetry -each evolution step is already gauge invariant and thus, although we obtain an approximate time evolution, the symmetry is exact. This was used in [16,17,26] for digital quantum simulation of lattice gauge theories (using a more "economical" ways of breaking the Hamiltonian to less pieces; the important thing is that intersecting links, involving common fermions, will not be included in the same Trotter step).

E. The Case of Quantum States: Gauging PEPS
We can also consider the states' point of view: one can try and construct the state from some local ingredients, that enable separate, individual actions of gauging transformations, and then, after gauging, project all the local pieces together in a way that forms a non-trivial, interacting state.
This can be done in the process we shall describe now. As we shall see, the states we construct in this way are nothing but PEPS -Projected Entangled Pair States [10] -but with a local symmetry, following the procedure of [14].
Step 1. At each site, construct a state |A (x) that will involve the physical matter degree of freedom, and some virtual degrees of freedom as well. The virtual degrees of freedom are associated with the edges of links starting or ending at the vertex, and out of them we will construct the virtual electric fields, E 0 (z), where z = r, u, l, d denotes the edge: right, up (outgoing) and left, down (ingoing). The state will be constructed in a way that a Gauss law will be satisfied by the physical charge and the virtual electric fields: i.e., we demand that it will be invariant under transformations generated by where Q is the charge associated with the physical matter degree of freedom at the vertex, which does not have to be fermionic. At each vertex we can define the virtual transformations and the physical one such that This allows us to write the symmetry condition of |A (x) in the conventional form used in the context of PEPS; at each vertex, the action with a group transformation on the physical degree of freedom is equivalent to acting on the virtual ones, Step 2. We take a product of the local matter states |A (x) with the gauge field states |0 on all the links starting at that vertex - Step 3. On each |A (x) we act with two gauging transformations U G , rotating (and entangling) the virtual degrees of freedom corresponding to links beginning at the vertex with respect to the gauge field degrees of freedom on that link: U G (x, 1) U G (x, 2). The local states obtained by that, will obey two symmetries, connecting the physical fields and the virtual fields on the outgoing links, These are generated by operators such as those we considered in the one-link case (30), i.e. connecting each physical electric field only to one rotated degree of freedom (the virtual one in this case) and thus gauging was possible. As a result, these states will also obey modified Gauss laws, in which the physical electric fields on the outgoing links replace the virtual ones, while the virtual fields remain for the ingoing links: These Gauss laws are already "half-way" to complete physical ones; we only need to contract the local states to one another, as done in the next step.
Step 4. We project the virtual ingredients into maximally entangled pairs connecting the edges of links |ω , satisfying the symmetry relations This will remove all the virtual degrees of freedom, after having used them for contracting the whole state, which is not a trivial product state anymore.
The final state Step 1

|0 x
Step 2 (58) will have the desired physical symmetry -local gauge invariance.
Steps 1 and 4 could easily be recognized as the construction of a PEPS: we project local states with some virtual symmetry to maximally entangled pair states. In fact, this is a PEPS as well -a gauged one. The lesson we learn from that, is that in order to obtain a PEPS description of a lattice gauge theory state, we can first write down a globally invariant PEPS, describing the matter degrees of freedom of the desired gauge theory, Step 1 (59) out of which the gauge invariant PEPS, with dynamical gauge fields, may be obtained in the process described above, that "interrupts" the usual PEPS construction in the form of steps 2,3. This is exactly what we discussed above: gauging, or minimally coupling, by a unitary transformation, must be done before we contract the local ingredients, that is, on a product state that allows us to act separately on different links. Only after gauging we project to a nontrivial physical state of the whole lattice, at step 4: |ψ is not gauged by directly acting with a unitary transformation on |ψ 0 ; rather, a modification is required before the system's ingredients are connected.
The nice "feature" that we get is that many symmetry properties of the original, globally invariant PEPS |ψ 0 "survive" the gauging procedure, such as spatial symmetries (translation, lattice rotation, lattice inversion etc.anything that might have existed prior to gauging) with the right modification for the gauge field transformation laws (e.g. translation invariance of |ψ 0 can be made charge conjugation symmetry of |ψ ). And, of course, the global symmetry becomes local. Therefore, it is enough to construct a state with a global symmetry, and the one with local gauge invariance is obtained immediately.
For the sake of completeness, and for readers not familiar with PEPS, let us show that the PEPS |ψ 0 ,|ψ indeed have the right symmetries. First, let us transform |ψ 0 with a global transformation and see what happens: -the transition from the first line to the second is thanks to the symmetry property of |A (x) ; from the third to the fourth, it is the symmetry property of the bond states ω |.
In the local case, As a final remark, note that this formulation can also be used for the construction of pure gauge states, with no physical matter: in that case, the states |A (x) will have no physical degrees of freedom, only virtual ones that enforce the Gauss law and connect the links.

Example
As an example [14], let us construct a very simple PEPS as follows. At each vertex, we use a physical Hilbert space spanned by states of the form {|p } J p=−J (eigenstates of Q, labeled by its eigenvalues) and virtual ones with {|v } j v=−j (similarly, eigenstates of the virtual electric fields E 0 ) where the integers j, J can be either finite or infinite. The symmetry is obtained for where |r, u, l, d is a product of four virtual basis states, {|v } -the Kronecker delta enforces the Gauss law (in general, the tensor elements A p ruld should be proportional to the right combination of Clebsch-Gordan coefficients that correspond to the right combination of representations [14]; in this simple U (1) case, the Clebsch-Gordan coefficients are just Kronecker deltas). The gauging transformation will simply be Note that as a result of the gauging that identifies the physical electric fields with the virtual ones, |E| ≤ j, and hence if we have a finite truncation for the virtual electric field, the physical electric field will be equally truncated as well.
Finally, the maximally entangled pair states on which we project will be

III. CONSTRUCTION OF GAUGE INVARIANT FERMIONIC PEPS (FPEPS)
While the above procedure for gauging a PEPS works for any PEPS with a global symmetry as |ψ 0 , we shall focus from now on on a more particular case, of gauged Gaussian fermionic PEPS. The free matter state with global symmetry we begin with is fermionic, in accordance with the usual matter formulations in gauge theories -and thus we will construct a fermionic PEPS (fPEPS) (for non fermionic lattice gauge PEPS, see [27] -with matter, and [28] -without).
Within the class of fermionic PEPS, the states we would like to construct and then gauge will be Gaussian fermionic PEPS [29]. Gaussian states are ground states of quadratic, free, non-interacting Hamiltonians, and thus gauging them is analogous to gauging a free matter theory through the conventional Hamiltonian (or Lagrangian) procedure. They are fully described in terms of their second moments (covariance matrix) -expectation values of fermionic bilinear, quadratic operators; all the other correlation functions are obtainable using the Wick theorem [30]. Therefore, they make sense as a starting point from both physical and computational grounds.
Once again, we will focus on the U (1) case in 2 + 1 dimensions, as it captures all the important qualitative features without mathematical complications that have to do with choosing some more complicated gauge group or higher dimensions. This derivation is discussed in detail in [13]. The SU (2) generalization may be found in [15], while a discussion for general gauge groups and dimensions is in [18].
We will hence begin with the construction of globally invariant fPEPS, and then turn to gauging them.

A. The Physical Setting
The physical ingredients we would like to include have already been discussed: a single fermionic mode, created by ψ † (x), is defined at each vertex. We will use a staggered fermionic formulation as before, but a slightly different one, which is obtained from the previous one by a particle-hole transformation of the odd sublattice. Previously, the absence of a fermion on an odd vertex represented the presence of an anti-particle. Here, the presence of a fermion will represent the presence of an anti-particle. This will allow us to define things in a more translationally invariant way on the level of states. The local fermionic charges are now and the Hamiltonian H f discussed above will transform (up to a constant) to the superconducting form, We wish to construct a state |ψ 0 with a global U (1) invariance; that is, |ψ 0 should be annihilated by or invariant under transformations generated by it, One can also demand further symmetries, such as translation and (lattice) rotation invariance, which we will discuss later on.

B. The Local Ingredients
Naively, we would first like to construct the states |A (x) , and their product ⊗ x |A (x) . However, a product state of fermions is not well-defined, since a tensor product factorization of a fermionic Fock space is not well defined. One could define some ordering but we would like to avoid that. Instead, at each vertex we will deal with an operator A (x), involving both the physical and virtual degrees of freedom, and construct the product state as where |Ω is some initial state of the system -including the Fock vacuum of the physical fermions. A product of fermionic operators is still not well defined without specifying some lattice ordering, unless all the A (x) operators commute. As each of them involves degrees of freedom defined on a different vertex, they will commute if and only if each A (x) has an even fermionic parity, which we shall thus demand. In order to do that, the virtual degrees of freedom must be represented by fermions as well, and we simply interpret |Ω as the total Fock vacuum, of both physical and virtual modes (otherwise, we must have an even local physical parity at each vertex, which does not make any physical sense; moreover, in our current example, with only one physical mode per vertex, even parity with non-fermionic virtual modes implies that the physical fermions cannot be excited). Now we can focus on a single vertex. Besides the physical mode ψ † we introduce eight virtual (auxiliary) fermionic modes, created by r † ± , u † ± , l † ± , d † ± , associated with the four legs attached to the vertex as before (right, up, left, down). Out of these, we construct the virtual electric fields, where z = r, u, l, d, and the virtual Gauss law operator -note that the staggering was incorporated into the definition of E 0 (z, x). The most general Gaussian state |A will be created from the Fock vacuum |Ω using operators of the form where α † i includes all the possible creation operators, either physical or virtual.
We wish to construct A (x) which are invariant under this virtual Gauss law operator, In order to do that, let us sort our modes to positive and negative ones, depending on the sign of the phase which is put on them by this transformation. The negative ones, a † i are the ones whose number operator appears in G 0 with a minus sign on an even vertex, and therefore while the positive ones, b † i are represented by number operators with a plus sign in front, giving rise to On an even vertex, the negative modes are a † It is easy to convince ourselves that the symmetry condition is satisfied if and only if the constructed state has the form Where T ij is a 5 × 4 matrix; the first row corresponds to the coupling virtual and physical modes, and the remaining square block, τ couples the virtual modes among themselves. A minimalistic choice is T ij (x) = T ij . Out of it we can construct a translationally invariant fPEPS (that will become charge conjugation invariant once we gauge): the staggering implies that when one translates by one lattice site, a † i ←→ b † i for the virtual modes (translation by one lattice site for staggered fermions is charge conjugation), and thus the square, virtual block τ has to be antisymmetric.
We obtain, for every x, that -the PEPS condition for global symmetry, as in the non-fermionic case.

C. Contracting the Globally Invariant fPEPS
Once again, to avoid the product of fermionic states, we use operators instead of states for the contraction. Since we wish our final state to be Gaussian, we will use Gaussian projectors for that. On a horizontal link, we define the (unnormalized) projector and on a vertical one, where in both cases Ω is the projector to the vacuum states of the modes on the link's edges. It is easy to see that for horizontal links while for vertical ones -this can be recognized as the usual projectors symmetry conditions for a PEPS, but different than the ones previously used, due to the staggering.
The fPEPS is then given by Step 4 x A (x) |Ω Step 1 where |Ω v is the virtual vacuum. The global transformations take the form e iΛQ = x even V (x) x odd V † (x) -but as this staggering has already been taken care of, e iΛQ |ψ 0 = |ψ 0 . To avoid confusing notations due to the staggering, we show this only graphically:

D. Gauging the PEPS
We go on to gauging the state |ψ . We introduce gauge field Hilbert spaces on the links, but how shall we relate them to the fermions? We simply have to modify the gauged PEPS we constructed previously to the fermionic case: Step 3 x A (x) Step 1

|Ω
(84) where U G ( ) is a gauging transformation that entangles the virtual fermion at the beginning of the link with the physical gauge field we introduce on it, i.e. it rotates the virtual fermion with respect to the gauge field. It is defined as We define what are the symmetry conditions satisfied by this object? It is easy to see that for an even vertex x, while for an odd one -the properties that allow us to gauge, just like in (30). From these two equations, and the symmetry condition of A (x), we obtain that for an even vertex, while for an odd one, This allows us to verify that |ψ is gauge invariant, i.e. that for an arbitrary x 0 , Once again we show it pictorially, for the two possible parities of the vertex around which the state is transformed. First, even, and finally odd, Another thing that has to be mentioned, is that in this method the physical electric field is truncated, and cannot take any integer eigenvalue, but rather only, in this case, 0, ±1. This is once again due to the fact that in the gauging process (the action of U G ) transfers the state of the virtual electric field to the physical one, and the virtual electric fields, constructed out of fermions in this particular manner (as the difference between two fermionic number operators) only take these values. In order to increase the truncation, one has to add more virtual fermions on each leg, such that higher eigenvalues of E 0 are obtained. This is called, in usual PEPS terms, increasing the bond dimension, and it also allows one, in general, to increase the number of parameters the PEPS depends on as well.

E. Other Symmetries
So far, we have constructed a fermionic Gaussian PEPS with a global U (1) symmetry, and made it local using our gauging procedure. What can we say about other symmetries obeyed by the state?
We have already commented how to encode translational invariance into the globally invariant state |ψ 0 , and that a translation by one lattice site, which is, for staggered fermions, charge conjugation, requires exchanging a † i ←→ b † i .
If we consider the staggering in the definition of U G ( ), we see that the transformation of the gauge field that leaves |ψ invariant is translation + inversion of sign, -which is exactly charge conjugation! Of course, it makes no sense to talk about simple translational invariance when we discuss staggered fermions, where the unit cell includes more than one vertex. If we consider translation by two sites, on the other hand, we have simple translational invariance, as expected for such a staggered system. Another symmetry that could be considered is rotation. First of all, for the globally invariant state without the gauge field. For that, one has to define the right transformation properties of all fermions under the rotation. This involves not only changing the coordinate respectively, but also requires more. As we discuss fermions, a complete rotation of the physical fermion must put a sign on it. Therefore, the basic π/2 rotation, realized by an operator U R , must give rise to with η 4 p = −1; Rx is the rotated coordinate. The virtual modes around each vertex must be exchanged by such a transformation, up to possible phases, keeping positive modes positive and negative modes negative. That is, we need to have, for the virtual components, with rotation matrices that contain only one nonzero element per row and column, and that has to be a phase. The phase must be chosen such that the product of projectors ω is invariant.
After gauging, one simply has to add the most obvious rotation rules for the gauge field. Once these two requirements (translation and rotation invariance) are taken care of, and we remove some redundancy having to do with some phase symmetries of the virtual fermions in which the physical fermions do not participate, we arrive at a final characterization [13], where t ≥ 0 and y, z ∈ C. t = 0 is of course meaningless in the global case, but in the local case it gives rise to pure gauge states.

F. Generalizations
One can generalize the gauging prescription for PEPS, and in particular Gaussian fPEPS, described above, in the following ways: 1. Different spatial dimensions: one simply has to introduce more similar virtual legs, and apply the gauging transformation only in the outgoing directions. A similar procedure holds for other geometries too.
2. Different gauge groups: similar procedure applies to further gauge groups (see [15] for an explicit SU (2) construction, and [18] for a general formulation). In general, for non-Abelian groups, the symmetry conditions become more complicated technically, as there is difference between left and right group transformations [12,22] (the electric field on the link differs from one edge to another -it has both left and right fields, whose difference manifests the fact that a non-Abelian gauge fields, unlike Abelian ones, carry charge -e.g. gluons vs. photons).
3. Larger bond dimensions: as explained before, one may add more virtual fermions in order to enlarge the dimension of the virtual Hilbert spaces and introduce more parameters on which the state depends. This, however, must be done carefully in a way that respects the symmetry (only complete irreducible multiplets can be included -and also multiple copies thereof). This increases the truncation of the virtual Hilbert spaces, and hence of the physical Hilbert space on the links as well (electric field truncation).

G. Analytical (Toy Model) Results
Which type of physics is manifested by "toy" states, with the minimal bond dimension (truncated electric fields), as described above?
One way to study that analytically is to put the system on a cylinder -with one compact and one open dimensionand define an object called the PEPS transfer matrix, similar to transfer matrix in statistical mechanics [31]. Roughly speaking, a two-point correlation function, for two operators L rows apart, depends on the L − 1 power of this transfer matrix; therefore, if the spectrum of the transfer matrix is gapped, correlation functions with respect to the PEPS show an exponential decay.
A conventional way to study such PEPS is therefore to compute the lowest eigenvalues of the transfer matrix for some different choices of the PEPS parameters, and see whether it is gapped or not -implying that the PEPS could be the ground state of a gapped Hamiltonian, or not. This allows one to draw a "phase diagram of states" in this parameter space -that is, a phase diagram for the set of Hamiltonians of which the set of PEPS obtained by varying the parameters are ground states (usually referred to as parent Hamiltonians). Parameter space surfaces on which the transfer matrix is gapless are suspected as phase boundaries.
In a case of a gauge invriant PEPS, such as ours, a two-point correlation function for two operators in two different rows, can be nonzero only if both operators are separately gauge invariant. Otherwise, a gauge field string must connect them. Therefore another type of a transfer matrix, with a U operator on it, is also physically meaningfulfor example, for the computation of the expectation value of mesonic operators, One could then draw a phase diagram in this method, and calculate the expectation values of some significant physical operators -order parameters -within the phases (gapped regions). Such studies, with Wilson loops [1], for U (1) states [13], have shown both confining and non-confining gapped phases for static charges, in the pure gauge case. Similar SU (2) studies [15] have shown a gapped deconfining phase as well as a gapless phase with perimeter-law decay of the area law -a hint of a possible Higgs phase [32].
Another way to detect possible phase transitions is to look for virtual symmetries: symmetries in which the physical degrees of freedom do not participate. In the SU (2) case [15], for example, such a particle-hole symmetry exists for the virtual fermions. This symmetry allows one to construct a parameter space transformation the that contexts states which are physically identical but constructed using different parameters. The fixed lines of this transformation are reasonably suspected as phase boundaries -and indeed, when those were compared to transfer matrix calculations, they gave rise to the same results.

H. Contracting Gauged Gaussian Fermionic PEPS: Combining Monte-Carlo with Tensor Network Methods
In many cases, tensor network states are used as ansatz states in variational calculations. While this has been done quite successfully for one dimensional systems, with DMRG (density matrix renormalization group) techniques [33], the higher dimensional generalizations are problematic as they do not scale well [10].
From the lattice gauge theory point of view, the widely used Monte-Carlo methods suffer from the sign problem [5] in several physically interesting scenarios (e.g. finite chemical potential). They also do not allow one to consider directly unitary time evolution, as they are formulated in a Euclidean spacetime.
Gaussian fermionic PEPS suggest a way to circumvent these problems. As we shall see, one may calculate, using Monte-Carlo methods, the expectation values of physical observables with respect to gauged Gaussian PEPS, exploiting the simple and efficient Gaussian formalism, based on Gaussian integration and matrix operations [30]. Such calculations do not depend on the dimension of the system, overcoming the first problem. They also do not encounter the sign problem, as they are carried out in the Hilbert space formulation (expectation values of physical operators) and not in the Wick rotated path integral formulation that encounters it.
For a general discussion see [18]; here, we shall briefly review the simple U (1) case.
We have already mentioned that the Hilbert space of a fermionic lattice gauge theory H, with no static charges, could be embedded in the product space where H g is the gauge field Hilbert space, covering all links, and H f is the fermionic Fock space, covering all vertices. The most general state in this Hilbert space may be expanded as where |Φ = ⊗ |φ ∈ H g is some gauge field configuration state on all the links, and DΦ = dφ ; |ψ (Φ) ∈ H f is a general fermionic state, that depends on the gauge field configuration as a parameter, which in general is not normalized. The state |ψ could be very general, not a PEPS and even not manifesting any symmetry; of course, we are interested in the case that it is gauge invariant. Then, |ψ (Φ) is a state of fermions coupled to an external, static U (1) field with configuration Φ. When |ψ is a gauged Gaussian fermionic PEPS, the states |ψ (Φ) are merely fermionic Gaussian states, and therefore admit a very simple description and analytical treatment using the Gaussian formalism [30]. This allows us to efficiently calculate the expectation values of physical operators for them.
Consider the (square) norm of the state. As the configuration state |Φ form a complete set, it can be simply expressed as -an integral over squares of norms of Gaussian states, ψ (Φ) |ψ (Φ) , whose computation is very simple. They are all positive definite, and thus the function is a probability density function in the space of gauge field configurations, and Z is the partition function associated with it. Now we are ready to compute expectation values. Consider first the Wilson loop [1], an oriented product of group element operators along some closed path C: (where is now an oriented link, that is, if its orientation is "backwards" one has to use U † on that link, such that the operator W (C) is gauge invariant). The configuration states are eigenstates of W (C), where now a minus sign is understood in front of the phase of a "backward" link. We therefore obtain that -which could be computed using Monte-Carlo. Let us consider another type of physical operator -a mesonic string, an open gauge field string enclosed by fermionic operators, e.g.
where now C is an open path connecting the endpoints x, y. In this case, we obtain that -this integral involves the expectation value of a quadratic fermionic operator with respect to a fermionic Gaussian -which is very simple to calculate in the Gaussian formalism using covariance matrix elements of |ψ (Φ) . Therefore, M (x, y, C) can be calculated too using Monte Carlo integration of quantities obtained from Gaussian calculations.
Similar results are obtained for other operators -such as ones involving electric fields [18]. In all cases, the quantum nature of the gauge field Hilbert space is handled before the integration takes place, leaving traces of the gauge field only in the form of integration variables. The quantum nature of the fermionic part is taken care of through the Gaussian formalism, thanks to which everything can be computed efficiently -as functions of the gauge field configurations. Eventually one is left with a DΦ integral over the gauge field configurations that could be calculated with Monte-Carlo. This opens the way to use such states, with higher bond dimensions, as varitational guesses for ground states of lattice gauge theory Hamiltonians, for example.

A. Solving Gauss Law for the Matter Field
We have argued that the Gauss law (42) is problematic to solve for the gauge field, since it is a differential equation with no unique solution, if the spatial dimension is more than one (and in the one dimensional case, what we get is nonlocal). However, note that if we see it as an equation for the matter, the setting is completely different: then, it is a very simple algebraic equation, in which the charge (or its density) -not a vector quantity -is already explicitly and uniquely given as a function of the divergence of the electric field.
Therefore, in cases where each local charge Q (x) eigenstate is obtained uniquely by a single configuration of the matter degree of freedom at x, one could indeed write down an "opposite" unitary transformation; that will take a gauge invariant state or Hamiltonian, with both matter and gauge fields as quantum degrees of freedom, and decouple the matter using the Gauss laws. This produces a state, or a Hamiltonian, without the symmetry, but also without the Gauss law constraints; with the same physical information, and effectively less degrees of freedom. It is, in fact, a controlled operation that changes the state of the matter at each vertex, controlled by the values of the electric field on the links around it. It is a well known procedure for Higgs matter fields, known as fixing the unitary gauge.

B. The Unitary Gauge of Higgs Theories
Let us consider a lattice formulation of the Higgs mechanism, similar to those of [32,52]. On each vertex, we have a complex scalar field Φ (x) which may be expanded in terms of two bosonic modes. We will write it in a polar form, In the usual quasi-classical treatment of the Higgs mechanism, R (x) is the Higgs field, and θ (x) is the Goldstone mode, which are independent degrees of freedom, satisfying The phase θ (x) is canonically conjugate to the local charge Q (x), having a non bounded integer spectrum, raised by e iθ(x) , and completely decopuled from the radial degree of freedom R (x). In [32,52], the Higgs field is frozen -R (x) = R 0 , and we will do it too, as it has no relevance for us. On the links, we have the usual U (1) gauge field Hilbert space introduced above (note that in this case, the link and matter Hilbert spaces are the same).
The interaction Hamiltonian between the matter and gauge fields takes the form As usual, in the Higgs mechanism discussion, we fix the gauge to the unitary gauge: perform a pure-gauge transformation (one that only acts on the gauge field degrees of freedom), in which the Goldstone modes are absorbed by the gauge field, making it massive. The transformation is How is it implemented quantum mechanically? It is a straightforward exercise to convince ourselves that it takes the form -where, by expressing it as a vertex-dependent transformation, we changed the point of view: instead of transforming the gauge field with respect to the matter configuration, we look at it now as a transformation that changes the matter field on each vertex, controlled by the electric fields on the links around it. It changes the matter charge at each vertex exactly by the amount of electric field divergence, and thanks to the Gauss law it means that it is reduced to zero. This is the controlled operation we wanted for decoupling the matter (while eliminating the gauge constraints)! If we act with it on a gauge invariant state in this Hilbert space, we will get a product of a trivial matter state with a gauge field superposition that still contains the same physical information, but without the symmetry: for an arbitrary gauge invariant state, with a Gauss law at each vertex, we have that without any Gauss law constraints. If we consider the case of a Z 2 gauge theory with such Higgs matter, where both the gauge field and the matter are described by spin half Hilbert spaces, U H , is nothing but a controlled not (CNOT) operation, as discussed in [27].

C. Eliminating Fermionic Matter
Next we wish to ask -is there an analogy for fermionic matter? The answer is: yes, but carefully: fermions come with their special statistics and its complications, and thus one has to be extra cautious, and act slightly differently than what we did above. Following [19], we will review (and demonstrate for U (1)) how to replace the fermionic degrees of freedom by spins (hard-core bosons).
We go back now to the U (1) lattice gauge theory with fermionic matter discussed in the very beginning, with interaction terms of the form In an (incomplete, as we shall see) analogy to the Higgs case, we express the fermionic mode operators as where c (x) is a fermionic Majorana mode, satisfying the Clifford algebra and the modes η † (x) have an on-site fermionic anticommutation relation, and off-site bosonic commutation relations, The operator η † (x) -almost analogous to the radial field R (x) -behaves like the creation operator of a hard-core boson; if we can remove the operators c (x) from the Hamiltonian, we can eventually represent η † (x) by spin raising operators, σ + (x). In this sense, the operators c (x) are analogous to e iθ(x) in the bosonic, Higgs case; indeed, one can think about the Majorana mode operator obeying c 2 (x) = 1 as a "Z 2 phase operator" having to do with the fermionic Z 2 symmetry. The analogy is not complete, though; in the Higgs case, the radial and angular components were independent, commuting on-site as well, but here, although commuting for different sites, on-site we have η † (x) , c (x) = 0: the "radial" and "angular" modes are related through a Z 2 equation. In fact, unlike in the Higgs case, the field η † (x) contains all the physical information, in the sense that its charge participates in the U (1) Gauss laws, while c (x) is only related to a Z 2 symmetry -the fermionic statistics. These two fields are related through the local Z 2 symmetry manifested by η † (x) , c (x) = 0. Therefore, in the fermionic case the "separation", which is not absolute, between "radius" and "angle" could be seen as separating physics and statistics -held together by the anticommutation relation.
In these terms, the interaction Hamiltonian has the form Ideally, we would like to have a transformation, for which and we are done. This, however, would be a problem: each link is transformed from the left by c (x) and from the right by c (x +ê i ). This is realized with local transformations of the form However, these transformations, that "rotate U (x, i) by a fermionic phase", do not have a fixed fermionic parity and violate the fermionic superselection locally. It can be shown that the product of all these, which is what we need, preserves fermionic parity globally, but as it is locally broken, one has to define some order for the product, that will eventually give rise to nonlocal results. What we do, instead, is to introduce on each = (x, i) link two Majorana modes, α i (x) associated with its beginning and β i (x +ê i ) with its end. These operators do not appear in the original Hamiltonian at all and hence we can fix the initial state of the respective auxiliary degrees of freedom as we wish. If the original Hilbert space, without these modes, is H, we simply embed it within H × Ω, where Ω is the one dimensional Hilbert space consisting of the Fock vacuum annihilated by the link annihilation operators as well as the vertex ones, χ (x), for which We define the local transformation It has an even fermionic parity, and thus preserves the fermionic superselection rule locally; thanks to that, and having different fermionic modes at each vertex, [U F (x) , U F (x )] = 0, and the transformation is well defined, since no ordering has to be specified for the product. We obtain that where ξ i are phases (signs) that depend on the electric fields on some links around the link's edges [19]. Since the initial state for the auxiliary link fermions was the vacuum, and one may check that [U F , α i (x) β i (x +ê i )] = 0, they remain so after the transformation, and could be replaced by a −i. The hard-core bosonic operators η † (x) can be effectively replaced by spin raising operators σ + (x) (in a process that involves local, string-less Jordan Wigner transformations -see [19]), and we obtain a transformed interaction Hamiltonian of the form Other terms in the Hamiltonian also transform in a way that leaves no fermions present [19]. Therefore we see that our lattice gauge theory with fermionic matter was mapped to one with hard core bosonic matter (that can be now removed straightforwardly using the Gauss laws). The lattice rotations are broken by the ξ i operators, whose role is to keep the right statistics and commutation relations [19]. In order to demonstrate it simply, let us look at the one dimensional case, where and one simply obtains As the "radial" and "angular" degrees of freedom were related to each other, the symmetry is not broken, and a matter field still exists; however, it is not fermionic. Unlike in the Higgs case, we did not break the continuous gauge symmetry completely and did not eliminate the matter: we only exploited the Z 2 subgroup of the gauge group, U (1) in our example, to eliminate the fermionic nature of the matter. In some cases, as the one we discuss, one may move on and eliminate the hard-core bosons using a simple, non-fermionic controlled operation as the one shown above for Higgs field; the difference will be, that now the bosons are hard-core, so projectors for the right values of electric field divergence must be introduced.
In fact, this can be done [19] for every lattice gauge theory whose gauge group includes a Z 2 normal subgroup -for example, Z 2N , U (N ) , SU (2N ): the important property is to have an operator E, for which the Z 2 relation e iπE U e −iπE = −U holds. This is the basis for our transformation. Groups with a Z 2 normal subgroup satisfy that. For other cases, e.g. SU (2N + 1) , one can introduce an auxiliary Z 2 gauge field -that is, extend the gauge symmetry to e.g. SU (2N + 1) × Z 2 , and use it for the elimination of fermions (and their replacement by hard core bosons) in a process as the one discussed above. It can be done in a way that preserves the physical properties of the original model (no dynamics is added for the auxiliary Z 2 gauge field). The relation to Z 2 does not look accidental: any "reasonable" fermionic theory without supersymmetry has a global Z 2 symmetry -parity superselection, which we have already discussed. When we have a Z 2 local symmetry, we can construct a transformation as above and eliminate the fermions locally since their information is "saved" by the gauge field. If there is no such field we can add it with some minimal coupling procedure and then do as above: so, in order to remove the fermionic statistics, one simply has to gauge the global Z 2 symmetry associated with it and make it local.