Projected entangled pair states for lattice gauge theories with dynamical fermions

TL;DR

The paper introduces a new method for simulating lattice gauge theories with fermions using projected entangled pair states, avoiding computational issues like the sign problem.

Contribution

The novel contribution is applying gauged Gaussian projected entangled pair states to lattice gauge theories with dynamical fermions.

Findings

The method shows agreement with exact diagonalization results for small systems.

The approach is computationally feasible for larger systems where exact solutions are not available.

This technique avoids the sign problem and can be extended to higher dimensions and other gauge groups.

Abstract

Lattice gauge theory is an important framework for studying gauge theories that arise in the Standard Model and condensed matter physics. Yet many systems (or regimes of those systems) are difficult to study using conventional techniques, such as action-based Monte Carlo sampling. In this paper, we demonstrate the use of gauged Gaussian projected entangled pair states as an ansatz for a lattice gauge theory involving dynamical physical matter. We study a \documentclass[12pt]{minimal}
\usepackage{amsmath}
\usepackage{wasysym}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{amsbsy}
\usepackage{mathrsfs}
\usepackage{upgreek}
\setlength{\oddsidemargin}{-69pt}
\begin{document}$${{\mathbb{Z}}}_{2}$$\end{document}Z2 gauge theory on a two dimensional lattice with a single flavor of fermionic matter on each lattice site. For small systems, our results show agreement with results computed by exactly diagonalizing the Hamiltonian, and demonstrate that the approach is computationally feasible for larger system sizes where exact results are unavailable. This is a further step on the road to studying higher dimensions and other gauge groups with manageable computational costs while avoiding the sign problem.

Simulating quantum field theories presents a host of computational challenges, including the large dimensionality of the problem and numerical obstacles such as the sign problem. The authors demonstrate a technique that is able to address these challenges, and show results for a lattice gauge theory which includes matter.

Figures

A diagram of the 2D lattice.Matter is shown on lattice sites in yellow and gauge fields on links in blue. The blue box shows the labeling convention for a plaquette.

A graphical representation of the construction of a gauged gaussian projected entangled pair state (GGPEPS) for a fixed gauge field configuration.We start on the left with the vacuum of the physical modes on each site (shown in yellow), create 4 virtual modes on each link around each site (turquoise), and couple them to each other and to the physical modes. We then couple the appropriate virtual modes to the gauge fields (blue). Finally, we project the virtual modes from neighboring sites onto a maximally entangled state, and trace out the virtual modes. The bottom row illustrates each operati

A graphical representation of the state which shows that an entanglement area law is automatically obeyed.The physical matter on the lattice sites is shown in yellow, the virtual modes are shown in turquoise, and the gauge fields are not shown. The entanglement between the region highlighted in purple and the rest of the lattice is proportional to the number of links through which the boundary cuts, since the only source of entanglement is the projection onto maximally entangled states of the virtual modes on the same link (enclosed by the turquoise dashed lines).

The various energy observables for the 2 × 2 ground states.The dots show results given byour ansatz, showing their close match to the exact diagonalization results (solid lines) for each term in the Hamiltonian. The interaction coupling *g*_I_ was fixed at *g*_I_ = 1.0, and the mass coupling *g*_M_ was fixed at *g*_M_ = 1.0. The horizontal axis is the electric coupling *g*_E_ = *λ*; the magnetic coupling is defined as *g*_B_ = 1/*λ*. The results shown here do not rely on Monte Carlo, but rather on exact contraction --- for a 2 × 2 system, it is possible to simply iterate over all $\documentcla

A parameter diagram for a 2 × 2 system with massless fermions (*g*_M_ = 0) coupled to gauge fields.The horizontal axis is the electric coupling *g*_E_ = *λ* (the magnetic coupling is given by 1/*λ*), and the vertical axis is the interaction coupling *g*_I_. The colors indicate the value of 1 × 1 Wilson loops for the ground state of the Hamiltonian with the specified couplings.

The ground state energies found for various lattice sizes and couplings.Moving from left to right, panels (a-b) show a 2 × 2 lattice (computed via exact contraction), panels (c-d) show a 4 × 4 lattice (computed with Monte Carlo), and panels (e-f) show a 6 × 6 lattice (computed with Monte Carlo). Where available, the results computed via exact diagonalization are shown as solid lines. The top row (panels a, c, e) is systems with *g*_M_ = 0 (massless fermions), while the bottom row (panels b, d, f) gives results for *g*_M_ = 1. The interaction coupling *g*_I_ is shown in the legend, while *λ* on

The expectation values of 2 × 2 Wilson loops for 4 × 4 and 6 × 6 systems.Though Wilson loops are not an order parameter for theories including matter, the sharp transition is noteworthy. The interaction coupling *g*_I_ is shown in the legend, while *λ* on the *x*-axis defines the electric and magnetic couplings (*g*_E_ = *λ* and *g*_B_ = 1/*λ*). The error bars shown in the legend (which are generally too small to be seen on the graph) are the errors due to Monte Carlo sampling.

Full text

Introduction

Gauge theories are central to the Standard Model of particle physics^1^, as well as condensed matter physics^2^. Yet such theories are notoriously hard to study. Discretizing space (or spacetime) onto a lattice—giving rise to a lattice gauge theory (LGT)—provides a standard framework for studying gauge theories, particularly using Monte Carlo methods^3–5^. However, doing so presents its own difficulties, prominent among them the sign problem, which arises when the probability intended for use in Monte Carlo does not form a valid probability distribution, i.e. takes on negative or complex values^6^.

We study a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{2}$$\end{document} lattice gauge theory with gauge fields on the links and fermionic matter on the lattice sites. Similar systems, including the pure \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{2}$$\end{document} gauge theory originally proposed by Wegner^7^ and models with Ising matter fields^8,9^, as well as fermionic matter fields, have been studied widely in several contexts. These include \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{2}$$\end{document} spin liquids^10,11^, unconventional quantum phases of matter^12–20^, and quantum simulation^21–27^. Here, we consider a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{2}$$\end{document} theory in two spatial dimensions with a single flavor of fermionic matter. This theory is a precursor to theories that suffer from a sign-problem and cannot be tackled with action-based Monte Carlo techniques.

Due to the above-mentioned computational problems, several methods have been proposed for the study of similar models in the last decade. One approach to the study of such systems is quantum simulation, with much recent progress and improvement as quantum hardware improves. Another is the use of tensor networks, and it is into this class which our work falls. Tensor networks are a class of ansatz states which are constructed by introducing extra “virtual” degrees of freedom, which are then traced out to produce a physical state^28,29^. They provide an important way of representing states based on entanglement properties, especially for states which obey an entanglement area law^30,31^. Algorithms based on tensor networks are an important tool, since their scaling is often polynomial in the size of the system under consideration^29,32–34^. Furthermore, tensor networks have been shown to capture the low energy states of many systems^30,32,33^, and they are useful for studying thermal states^35^, as well as time evolution^36,37^.

The properties of tensor network states in one dimension, known as matrix product states (MPS), are relatively well-understood^29^. In higher dimensions however, such as when working with projected entangled pair states (PEPS), which are the higher dimensional generalization of MPS, difficulties arise with both the theoretical and numerical aspects^38^. These issues largely have to do with the computational scaling of contracting the tensor networks and computing observables^29,39,40^.

Lattice gauge theories with tensor networks have been studied using several approaches. In a single space dimension, tensor network methods have been shown to overcome the sign problem and simulate time evolution well; details may be found in the review papers^41,42^ and references therein. Some numerical work has also been done in higher dimensions, such as the pure-gauge works^43,44^, as well as works using infinite PEPS (iPEPS)^45,46^ and tree tensor networks (TTN)^47–51^.

While these methods are very powerful numerically, they are not especially tailored to physical gauge symmetries. PEPS, with their ability to describe global symmetries using a virtual local symmetry^29^ suggest another path for studying lattice gauge theories in higher dimensions using tensor network states, building upon the concept of gauging, introduced in refs. ^52^, ^53^, through which globally invariant PEPS are made physically gauge invariant locally. We focus on the gauging procedure introduced in the latter reference, which was proven to be general enough for lattice gauge theory PEPS^54,55^, and in particular for constructing gauged Gaussian projected entangled pair states. Note that another symmetry-tailored method is that of tensor renormalization group (see, e.g., the review^56^), but it focuses on path integral methods while we aim here at the Hamiltonian formalism with tensor network states.

Gauged Gaussian projected entangled pair states (GGPEPS) were introduced as a class of ansatz states designed for lattice gauge theories^57,58^. They build on the work of^59^ which defined (ungauged) Gaussian PEPS. GGPEPS provide a convenient way to ensure that the local constraints following from gauge invariance are met. As PEPS, they satisfy an entanglement area law—the entanglement between two subsystems is proportional to the area (or lower/higher dimensional analogues) of the boundary between them. This is valuable, as the ground states of many LGTs of interest are known (in 1D) or conjectured (in higher dimensions) to obey such a law^31^. GGPEPS also allow for the efficient computation of observables, including standard observables such as Wilson loops or mesonic operators as well as gradients of operators, which is important for calculating the gradients of the Hamiltonian for a ground state search via a minimization procedure^60^.

In this paper, we present results using gauged Gaussian PEPS for a theory with dynamical fermionic matter. In reference^61^ the theory surrounding these states was developed in a general, abstract, and rigorous fashion, and in refs. ^62^^,^^63^, these states were used to investigate the ground states of pure-gauge theories. This paper improves upon those by demonstrating numerical results with dynamic fermionic matter. While the particular theory we consider does not suffer from the sign problem, it serves as a further demonstration that this approach offers a viable method for studying such theories as we work towards models that do suffer from the sign problem.

In the next section we describe the physical system we use for illustration—a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{2}$$\end{document} lattice gauge theory with dynamical fermionic matter. In the following section, we present our ansatz state and algorithm for finding the ground state. In the Results section we present our results, including arguments that the results capture the physics of the desired state. We conclude with some reflections and plans for future work in the Discussion section.

Methods

Physical system

A lattice gauge theory, such as the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{N}$$\end{document} one we study here, has matter (fermionic, in our case) on the lattice sites and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{N}$$\end{document} gauge fields on the links. We consider a two dimensional lattice of size Lx by Ly with periodic boundary conditions, and describe each lattice site by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{x}}=(x,y)\in {{\mathbb{Z}}}^{2}$$\end{document} and each link by ℓ = (x, k), where k ∈ {1, 2} denotes the direction of the link attached to the site x. This is illustrated in Fig. 1. We label unit vectors in each direction by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widehat{{\bf{e}}}}_{1}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widehat{{\bf{e}}}}_{2}$$\end{document} . It is occasionally convenient to refer to all the links around a given site; we therefore allow k ∈ {1, 2, 3, 4} and identify link ℓ = (x, 3) with the link \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$({\bf{x}}-{\widehat{{\bf{e}}}}_{1},1)$$\end{document} for horizontal links, and ℓ = (x, 4) with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$({\bf{x}}-{\widehat{{\bf{e}}}}_{2},2)$$\end{document} for vertical links. It will be necessary to treat sites on the even and odd sublattices differently; we therefore define (−1)^x^ = (−1)^x+y^ where x, y are the coordinates of the site x, so that (−1)^x^ = 1 on the even sublattice and (−1)^x^ = − 1 on the odd sublattice.Fig. 1A diagram of the 2D lattice.Matter is shown on lattice sites in yellow and gauge fields on links in blue. The blue box shows the labeling convention for a plaquette.

The gauge fields are represented by group elements \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g\in {{\mathbb{Z}}}_{N}$$\end{document} , which we use to construct an N-dimensional Hilbert space on each link. On each such Hilbert space we choose two unitary operators, P and Q, which satisfy

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${P}^{N} = {Q}^{N}={\mathbb{1}}\\ {P}^{\dagger }P = {Q}^{\dagger }Q={\mathbb{1}}\\ PQ{P}^{\dagger } = {e}^{i\delta }Q,$$\end{document}

where δ = 2π/N. We denote the eigenstates of P by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{\left|p\right\rangle \}$$\end{document} , which have eigenvalue e^ipδ^, and the eigenstates of Q by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{\left|q\right\rangle \}$$\end{document} which have eigenvalue e^iqδ^. Q is then a raising operator on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{\left|p\right\rangle \}$$\end{document} while P is a lowering operator on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left|q\right\rangle$$\end{document} , where the basis states are periodic, wrapping back to zero at p = N or q = N ^64^.

We will focus here on the case of a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{2}$$\end{document} theory to demonstrate the method. In this case, these operators can be represented by Pauli matrices, for example P = σz = P^†^ and Q = σx = Q^†^.

We consider a single flavor of fermionic matter, without spin or color; to reduce the doubling problem^65^ we stagger^66^ — placing the two-components of the matter/anti-matter spinor on even/odd lattice sites respectively. A fermion on site x is created by the operator ψ^†^(x), which obeys the anti-commutation relation for fermions, {ψ^†^(x), ψ(y)} = δ(x, y), where δ is the Kronecker delta.

The system is described by a Hamiltonian composed of terms which act on the gauge fields alone, the matter alone, and on both together^4^. For a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{N}$$\end{document} gauge theory, the terms which depend solely on the gauge fields are given by ref. ^64^ as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${H}_{{\rm{E}}} = \mathop{\sum }\limits_{\ell }(2-({P}_{\ell }+{P}_{\ell }^{\dagger }))\\ {H}_{{\rm{B}}} = \mathop{\sum }\limits_{p}2(1-{Q}_{1}{Q}_{2}{Q}_{3}^{\dagger }{Q}_{4}^{\dagger }),$$\end{document}

and for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{2}$$\end{document} specifically, this becomes

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${H}_{{\rm{E}}} =\mathop{\sum }\limits_{\ell }2(1-{P}_{\ell })\\ {H}_{{\rm{B}}} =\mathop{\sum }\limits_{p}2(1-{Q}_{1}{Q}_{2}{Q}_{3}{Q}_{4}),$$\end{document}

since P and Q are Hermitian in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{2}$$\end{document} theory. The offsets ensure that these terms are positive; they do not affect the physics. The subscripts E and B stand for “electric” and “magnetic”, where these terms are borrowed from the U(1) theory of electromagnetism. In the U(1) theory, the analogous terms have that meaning (and in the large N limit, compact QED, a lattice U(1) model^67^, is obtained from the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{N}$$\end{document} Hamiltonian^64^). As in that case, the electric term acts locally on the gauge fields, while the magnetic term is related to loops. The first term is a sum over all links, while the second is over all plaquettes (we have labelled the Q operators by their position in the plaquette rather than with links labelled by site and direction, as shown in Fig. 1).

Next, we add a term acting solely on the matter degrees of freedom

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${H}_{{\rm{M}}} = \frac{{L}_{x}{L}_{y}}{2}+\mathop{\sum }\limits_{{\bf{x}}}{(-1)}^{{\bf{x}}}{\psi }^{\dagger }({\bf{x}})\psi ({\bf{x}})\\ = \mathop{\sum }\limits_{{\bf{x}}}\left(\frac{1}{2}+{(-1)}^{{\bf{x}}}{\psi }^{\dagger }({\bf{x}})\psi ({\bf{x}})\right),$$\end{document}

where again the offset bounds the term below by zero. These are staggered fermions^66^, chosen for simplicity, where matter is placed on the even sites and anti-matter on the odd sites.

Finally we introduce an interaction term, which acts on each pair of neighboring sites together with the gauge field between them,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${H}_{{\rm{I}}} = \mathop{\sum }\limits_{{\bf{x}}}\left({\eta }^{2}{\psi }^{\dagger }({\bf{x}})U({\bf{x}},1)\psi ({\bf{x}}+{\widehat{{\bf{e}}}}_{1})+h.c.\right)\\ -\mathop{\sum }\limits_{{\bf{x}}}{(-1)}^{{\bf{x}}}\left({\psi }^{\dagger }({\bf{x}})U({\bf{x}},2)\psi ({\bf{x}}+{\widehat{{\bf{e}}}}_{2})+h.c.\right).$$\end{document}

The phases η = e^i**π/4^ and (−1)^x^ relate to rotations, as discussed below. The gauging operator U which acts on the gauge fields is defined for a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{N}$$\end{document} gauge theory as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$U({\bf{x}},k)=\mathop{\sum }\limits_{0\le q < N}{e}^{i\delta q}\left|q\right\rangle \left\langle q\right|,$$\end{document}

which gives that U(x, k) = Q(x, k) = σx(x, k) for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{2}$$\end{document} . For a general group other than \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{N}$$\end{document} , the system is described by a Hamiltonian given in refs. ^4,68^.

Unlike the other terms, the interaction energy is not bounded below by zero. Furthermore, it is not trivial to find it’s minimal value, which in general is not proportional to the lattice size. For a 2 × 2 lattice, the lowest possible value of this term can be calculated analytically as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-4\sqrt{2}$$\end{document} . For larger systems, the minimal value can be found using band structure considerations, as detailed in Supplementary Information C: Free Fermion Limit.

The full Hamiltonian is thus

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H={g}_{{\rm{E}}}{H}_{{\rm{E}}}+{g}_{{\rm{B}}}{H}_{{\rm{B}}}+{g}_{{\rm{I}}}{H}_{{\rm{I}}}+{g}_{{\rm{M}}}{H}_{{\rm{M}}},$$\end{document}

where gK are the couplings of the individual terms. Since there is only one flavor of matter, we do not add a chemical potential term. We denote the electric coupling gE by λ, and set the magnetic coupling to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{1}{\lambda }$$\end{document} , in order to match the relationship that arises through the discretization procedure of a continuous quantum field theory (QFT) onto a square lattice. Without consideration of the origins of the Hamiltonian in QFT, the couplings could be taken to vary independently.

The Hamiltonian and the states in which we are interested are invariant under several symmetries. The most central is local gauge invariance under the operators

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V({\bf{x}})={(-1)}^{{\bf{x}}}{P}_{{\bf{x}},1}{P}_{{\bf{x}},2}{P}_{{\bf{x}},3}{P}_{{\bf{x}},4}{e}^{i\pi {\psi }^{\dagger }({\bf{x}})\psi ({\bf{x}})},$$\end{document}

i.e. [H, V(x)] = 0 for all lattice sites x. Note that we have here used the fact that P = P^†^ for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{2}$$\end{document} . This symmetry divides the Hilbert space into sectors which can be labeled by the eigenvalues of the operators V(x). We study the sector satisfying

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V({\bf{x}})\left|\psi \right\rangle =\left|\psi \right\rangle \,\forall \,{\bf{x}},$$\end{document}

i.e. the sector with no static charges anywhere on the lattice.

Since there are no static charges, the system is invariant under translations by two sites — it is not invariant under translations of one site due to the staggering, which leads to the factor of (−1)^x^ in Equation (8).

The Hamiltonian is also invariant under rotations; since we work on a square lattice, the allowed rotations are those which rotate the plane by multiples of π/2. Since fermions pick up a minus sign upon a full rotation, we must account for this in the prescription for the rotation of fermionic operators,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\psi }^{\dagger }({\bf{x}})\to \left\{\begin{array}{l}\eta {\psi }^{\dagger }(\Lambda {\bf{x}})\,{\bf{x}}\,\,{\rm{even}}\,\\ \bar{\eta }{\psi }^{\dagger }(\Lambda {\bf{x}})\,{\bf{x}}\,\,{\rm{odd}}\,\end{array}\right.$$\end{document}

where Λx = Λ(x, y) = ( − y, x) implements the rotation, and η = e^i**π/4^. The operators on the links do not pick up a phase, and transform as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q({\bf{x}},k)\to \left\{\begin{array}{ll}Q(\Lambda {\bf{x}},2) & k=1\\ Q(\Lambda {\bf{x}}-{\widehat{{\bf{e}}}}_{1},1) & k=2.\end{array}\right.$$\end{document}

Only with rotations defined in this way is the Hamiltonian invariant under the allowed rotations.

Finally, the Hamiltonian is invariant under a global U(1) symmetry, ψ^†^(x) → e^i**θ^ψ^†^(x) for all sites. This corresponds to a total number conservation, i.e. the symmetry is generated by the conserved quantity

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=\mathop{\sum }\limits_{{\bf{x}}}{\psi }^{\dagger }({\bf{x}})\psi ({\bf{x}}).$$\end{document}

Since [H, N] = 0, this operator also divides the Hilbert space into sectors labelled by its eigenvalues; we focus on the “half-filling” case, in which there are equal numbers of matter and anti-matter particles. This sector is that of the Dirac sea, in which all even sites are empty while odd sites are full, and is given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\left|\psi \right\rangle =\frac{{L}_{x}{L}_{y}}{2}\left|\psi \right\rangle .$$\end{document}

We are interested in finding the ground state energy of the Hamiltonian given in Equation (7) and describing its physics. To do so, we build an ansatz state and preform a variational Monte Carlo minimization over its variational parameters, as we describe in the next section. Once a particular state of interest—in our case, the ground state—has been found, one may also be interested in observables other than the energy. Particular observables of interest include Wilson and Polyakov loops as well as mesonic operators. The closed loop operators take the form of products of link operators U along the path determined by the loop; mesonic operators are similar, but with matter operators ψ^†^ and ψ at the ends of the string. In both cases, the loop or string is a directed path, and one must consider the conjugate loop operators on those links which are traversed in the negative directions. In the case of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{2}$$\end{document} theory studied here, this slight complication can be ignored, since all loop operators are Hermitian. In the Results section we show some examples of Wilson loops calculations on the ground state.

Ansatz and algorithm

Particle hole transformation

Before constructing the ansatz state, it is useful to perform a transformation that will simplify the construction, and make the resulting system and ansatz state translationally invariant under translations by a single site. The transformation swaps the creation and annihilation operators on the odd sublattice, which allows us to treat the anti-matter on the odd sites in exactly the same way as the matter on the even sites. We therefore call it a “particle-hole transformation.” It is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\psi }^{\dagger }({\bf{x}}) \to \left\{\begin{array}{ll}{\widetilde{\psi }}^{\dagger }({\bf{x}}) & {\bf{x}}\,\,{\rm{even}}\\ \widetilde{\psi }({\bf{x}}) & {\bf{x}}\,\,{\rm{odd}}\end{array}\right.\\ \psi ({\bf{x}}) \to \left\{\begin{array}{ll}\widetilde{\psi }({\bf{x}}) & {\bf{x}}\,\,{\rm{even}}\\ {\widetilde{\psi }}^{\dagger }({\bf{x}}) & {\bf{x}}\,\,{\rm{odd}}\,.\end{array}\right.$$\end{document}

Note that this mathematical transformation does not change the underlying physics.

After the transformation, the system is invariant under translations by a single site. The symmetry under rotations remains, but with the modification of Equations (10) and (11) to

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{array}{rcl}{\widetilde{\psi }}^{\dagger }({\bf{x}}) & \to & \quad \eta {\widetilde{\psi }}^{\dagger }(\Lambda {\bf{x}}),\hfill\\ Q({\bf{x}},k) & \to & \left\{\begin{array}{ll}Q(\Lambda {\bf{x}},2)\hfill & k=1\\ -Q(\Lambda {\bf{x}}-{\widehat{{\bf{e}}}}_{1},1) & k=2.\end{array}\right.\end{array}$$\end{document}

Since U(x, k) = Q(x, k), this also gives the transformation rules for U(x, k). The global U(1) symmetry of total particle number conservation is also transformed, into

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\widetilde{\psi }}^{\dagger }({\bf{x}}){\to }_{U(1)}\left\{\begin{array}{ll}{e}^{i\theta }{\widetilde{\psi }}^{\dagger }({\bf{x}}) & {\bf{x}}\,\,{\rm{even}}\\ {e}^{-i\theta }{\widetilde{\psi }}^{\dagger }({\bf{x}}) & {\bf{x}}\,\,{\rm{odd}}\,,\end{array}\right.$$\end{document}

while gauge invariance is now defined as invariance under

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V({\bf{x}})={P}_{{\bf{x}},1}{P}_{{\bf{x}},2}{P}_{{\bf{x}},3}^{\dagger }{P}_{{\bf{x}},4}^{\dagger }{e}^{i\pi {\widetilde{\psi }}^{\dagger }({\bf{x}})\widetilde{\psi }({\bf{x}})}.$$\end{document}

Finally, we must rewrite the Hamiltonian in terms of the transformed fermionic operators. The operators which act only on the gauge fields remain unchanged, but the mass and interaction terms become

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${H}_{{\rm{M}}} = \mathop{\sum }\limits_{{\bf{x}}}{\widetilde{\psi }}^{\dagger }({\bf{x}})\widetilde{\psi }({\bf{x}}),\hfill\\ {H}_{{\rm{I}}} = \mathop{\sum }\limits_{{\bf{x}}}\left[{\eta }^{2}{\widetilde{\psi }}^{\dagger }({\bf{x}})Q({\bf{x}},1){\widetilde{\psi }}^{\dagger }({\bf{x}}+{\widehat{{\bf{e}}}}_{1})\right. \left.+{\bar{\eta }}^{2}\widetilde{\psi }({\bf{x}}+{\widehat{{\bf{e}}}}_{1})Q({\bf{x}},1)\widetilde{\psi }({\bf{x}})\right]\\ -\left[{\widetilde{\psi }}^{\dagger }({\bf{x}})Q({\bf{x}},2){\widetilde{\psi }}^{\dagger }({\bf{x}}+{\widehat{{\bf{e}}}}_{2})\right. \left.+\widetilde{\psi }({\bf{x}}+{\widehat{{\bf{e}}}}_{2})Q({\bf{x}},2)\widetilde{\psi }({\bf{x}})\right].$$\end{document}

In the remainder of this paper, we drop the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{\cdot }$$\end{document} on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\widetilde{\psi }$$\end{document} .

Ansatz

A general state in the Hilbert space of the full system — matter and gauge fields together — is a ray that can be written

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left|\psi \right\rangle =\mathop{\sum }\limits_{F,{\mathcal{G}}}\varphi (F,{\mathcal{G}})\left|F\right\rangle \left|{\mathcal{G}}\right\rangle ,$$\end{document}

i.e. as a sum over basis states of the fermionic matter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left|F\right\rangle$$\end{document} and gauge fields \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left|{\mathcal{G}}\right\rangle$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi (F,{\mathcal{G}})$$\end{document} is a complex amplitude. However, in general, such a state will not be gauge invariant — to ensure gauge invariance one must impose conditions on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi (F,{\mathcal{G}})$$\end{document} which rule out non-gauge-invariant states. Defining the gauge field configuration \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left|{\mathcal{G}}\right\rangle$$\end{document} as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left|{\mathcal{G}}\right\rangle {=\bigotimes }_{{\bf{x}},k}\left|g\left({\bf{x}},k\right.\right\rangle ,$$\end{document}

one can always write a state in the form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left|\psi \right\rangle =\mathop{\sum }\limits_{{\mathcal{G}}}\,{\psi }_{I}({\mathcal{G}})\left|{\psi }_{II}({\mathcal{G}})\right\rangle \left|{\mathcal{G}}\right\rangle ,$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left|{\psi }_{II}({\mathcal{G}})\right\rangle$$\end{document} is a state of the matter, which now depends on the gauge field configuration in a way that ensures gauge invariance, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\psi }_{I}({\mathcal{G}})$$\end{document} is a complex amplitude. The sum is over all possible gauge field configurations on all the links, which we express in the basis of eigenstates of U(x, k) = Q(x, k), denoted the magnetic basis (though Equation (21) remains valid in any basis).

We wish to construct an ansatz for both \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\psi }_{I}({\mathcal{G}})$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left|{\psi }_{II}({\mathcal{G}})\right\rangle$$\end{document} which satisfies the desiderata already mentioned: it should be computationally tractable to evaluate observables, gauge invariant, and result in a final state which obeys an entanglement area law. The GGPEPS ansatz meets all of these conditions.

We construct \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\psi }_{I}({\mathcal{G}})$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left|{\psi }_{II}({\mathcal{G}})\right\rangle$$\end{document} as fermionic Gaussian PEPS^69^ in very similar ways (in the language of^61^, they are simply different layers of the ansatz), but with the inclusion of physical matter for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left|{\psi }_{II}({\mathcal{G}})\right\rangle$$\end{document} . A graphical representation of the construction is shown in Fig. 2.Fig. 2A graphical representation of the construction of a gauged gaussian projected entangled pair state (GGPEPS) for a fixed gauge field configuration.We start on the left with the vacuum of the physical modes on each site (shown in yellow), create 4 virtual modes on each link around each site (turquoise), and couple them to each other and to the physical modes. We then couple the appropriate virtual modes to the gauge fields (blue). Finally, we project the virtual modes from neighboring sites onto a maximally entangled state, and trace out the virtual modes. The bottom row illustrates each operation locally; it is in fact applied to the entire lattice at each stage.

As in many tensor-network approaches, we begin by introducing virtual modes — for each site, in each direction, we add new degrees of freedom. We label these modes by the link direction with which they are associated k ∈ {1, 2, 3, 4} as well as an index μ to indicate the “copy”, since we add multiple virtual modes per link per site. In our case, we add four modes per link per site (for reasons discussed in Supplementary Information A: Ansatz Details), giving a total of 16 virtual modes per site created by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\alpha }_{k\mu }^{\dagger }({\bf{x}})$$\end{document} . For convenience, we will often denote the modes with a single index that wraps k and μ, and will include the physical mode ψ^†^(x) among them.

The new virtual modes are coupled to each other and to the physical matter associated with the appropriate site through an operator of the Gaussian form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A({\bf{x}})=\exp \left({{\mathcal{T}}}_{ij}({\bf{x}}){\alpha }_{i}^{\dagger }({\bf{x}}){\alpha }_{j}^{\dagger }({\bf{x}})\right),$$\end{document}

where i, j run over all the physical and virtual modes, which we denote by α. Summation is implied over i, j. This operator acts on the product of the physical and virtual Fock vacua \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left|{\Omega }_{{\rm{p}}}\right\rangle ,\left|{\Omega }_{{\rm{v}}}\right\rangle$$\end{document} . Since there is one physical mode per site, and we have added 4(2d) = 16 (where d = 2 is the dimension of the lattice) virtual modes, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal{T}}$$\end{document} is a 17 × 17 matrix. In order to ensure the final state obeys all the symmetries of the system discussed above, we must impose constraints on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal{T}}$$\end{document} — the details of how to do so are in Supplementary Information A: Ansatz Details. Of the four virtual modes per site per link, two will be used in the construction of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\psi }_{I}({\mathcal{G}})$$\end{document} and two in the construction of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left|{\psi }_{II}({\mathcal{G}})\right\rangle$$\end{document} . These sets of modes will not couple to each other, and the physical mode will only couple to the virtual modes used in the construction of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left|{\psi }_{II}({\mathcal{G}})\right\rangle$$\end{document} , so A(x) can be decomposed as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A({\bf{x}})={A}_{I}({\bf{x}}){A}_{II}({\bf{x}}),$$\end{document}

where each AJ(x) has the form given in Equation (22), involving only some of the modes — for AI(x), two copies of virtual modes per link (for a total of 8); for AI**I(x), two other copies of virtual modes per link plus the single physical mode on the site x (for a total of 9).

Next we introduce a gauging operator, following the procedure of^53^ which was proven to be general in^54,55^. This operator couples the gauge field on each link to the virtual modes just introduced. The coupling guarantees that the state will be gauge invariant. It is sufficient to only couple the gauge field on a link to the virtual modes associated with one of the neighboring sides; we couple the field on each link to the virtual modes associated with the site x when the link is represented as ℓ = (x, k) with k ∈ {1, 2}, i.e. the site to the left or below the given link. For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{2}$$\end{document} , the gauging operator can be written as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{U}}}^{{\mathcal{G}}}({\bf{x}},k)={\mathbb{1}}\otimes {{\mathbb{P}}}_{{\rm{even}}}+Q({\bf{x}},k)\otimes {{\mathbb{P}}}_{{\rm{odd}}}$$\end{document}

where here αi ranges over the virtual modes of link (x, k) which arise from x (but not those which arise from site \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bf{x}}\pm {\widehat{{\bf{e}}}}_{k}$$\end{document} ). In writing Equation (24), we have made use of Equation (6), which gives that U(x, k) = Q(x, k). \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{P}}}_{{\rm{even}}}$$\end{document} is a projection operator onto the state with an even number of excited virtual modes; \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{P}}}_{{\rm{odd}}}$$\end{document} is defined analogously. This gauging operator can also be decomposed into operators which act only on the virtual modes used in the construction of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\psi }_{I}({\mathcal{G}})$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left|{\psi }_{II}({\mathcal{G}})\right\rangle$$\end{document} , so that

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{U}}}^{{\mathcal{G}}}({\bf{x}},k)={{\mathcal{U}}}_{I}^{{\mathcal{G}}}({\bf{x}},k){{\mathcal{U}}}_{II}^{{\mathcal{G}}}({\bf{x}},k),$$\end{document}

where each of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{U}}}_{J}^{{\mathcal{G}}}({\bf{x}},k)$$\end{document} are defined as in Equation (24) over the appropriate virtual modes. Note that both \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathcal{U}}}_{J}^{{\mathcal{G}}}({\bf{x}},k)$$\end{document} operate on the same — physical — gauge field.

The final step in the construction of the ansatz state is the projection of the virtual modes from neighboring sites (on the same link) onto a maximally entangled state, and tracing out the virtual modes to leave a state with only physical degrees of freedom (the physical matter and the gauge fields). We define

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w({\bf{x}},k)=\exp \left({W}_{ij}^{(k)}{\alpha }_{i}^{\dagger }({\bf{x}}){\alpha }_{j}^{\dagger }({\bf{x}}+{\widehat{{\bf{e}}}}_{k})\right)$$\end{document}

with summation implied over i, j. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${W}_{ij}^{(k)}$$\end{document} determines the coupling between the virtual modes. It is convenient to define notation for the virtual modes such that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${r}_{\mu }^{\dagger },{u}_{\mu }^{\dagger },{l}_{\mu }^{\dagger },{d}_{\mu }^{\dagger }$$\end{document} as those which are right, up, left, down of the site, as shown in Fig. 2. Thus

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{array}{rcl}{r}_{\mu }^{\dagger }({\bf{x}}) = {\alpha }_{1\mu }^{\dagger }({\bf{x}})\,\qquad & {u}_{\mu }^{\dagger }({\bf{x}})={\alpha }_{2\mu }^{\dagger }({\bf{x}})\\ {l}_{\mu }^{\dagger }({\bf{x}}) = {\alpha }_{3\mu }^{\dagger }({\bf{x}})\,\qquad & {d}_{\mu }^{\dagger }({\bf{x}})={\alpha }_{4\mu }^{\dagger }({\bf{x}}).\\ \end{array}$$\end{document}

For the particular ansatz used in the results below, we choose W^(k)^ to give

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${w}_{I}({\bf{x}},1) = \exp \left({r}_{1}({\bf{x}}){l}_{2}({\bf{x}}+{\widehat{{\bf{e}}}}_{1})+{r}_{2}({\bf{x}}){l}_{1}({\bf{x}}+{\widehat{{\bf{e}}}}_{1})\right)\hfill\\ {w}_{I}({\bf{x}},2) = \exp \left({\eta }^{2}{u}_{1}({\bf{x}}){d}_{2}({\bf{x}}+{\widehat{{\bf{e}}}}_{2})+{\eta }^{2}{u}_{2}({\bf{x}}){d}_{1}({\bf{x}}+{\widehat{{\bf{e}}}}_{2})\right)$$\end{document}

for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\psi }_{I}({\mathcal{G}})$$\end{document} , and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${w}_{II}({\bf{x}},1) = \exp \left({r}_{3}({\bf{x}}){l}_{3}({\bf{x}}+{\widehat{{\bf{e}}}}_{1})+{r}_{4}({\bf{x}}){l}_{4}({\bf{x}}+{\widehat{{\bf{e}}}}_{1})\right)\hfill\\ {w}_{II}({\bf{x}},2) = \exp \left({\eta }^{2}{u}_{3}({\bf{x}}){d}_{3}({\bf{x}}+{\widehat{{\bf{e}}}}_{2})+{\eta }^{2}{u}_{4}({\bf{x}}){d}_{4}({\bf{x}}+{\widehat{{\bf{e}}}}_{2})\right)$$\end{document}

for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left|{\psi }_{II}({\mathcal{G}})\right\rangle$$\end{document} , which together give

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$w({\bf{x}},k)={w}_{I}({\bf{x}},k){w}_{II}({\bf{x}},k).$$\end{document}

The symmetries under which the projecting operators are invariant are discussed in Supplementary Information A: Ansatz Details; see also ref. ^61^.

With all of these ingredients we define

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\psi }_{I}({\mathcal{G}})=\left\langle {\Omega }_{I}\right|\mathop{\prod }\limits_{{\bf{x}},k}{w}_{I}({\bf{x}},k)\mathop{\prod }\limits_{{\bf{x}},k}{{\mathcal{U}}}_{I}^{{\mathcal{G}}}({\bf{x}},k)\mathop{\prod }\limits_{{\bf{x}}}{A}_{I}({\bf{x}})\left|{\Omega }_{I}\right\rangle$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{array}{rcl}\left|{\psi }_{II}({\mathcal{G}})\right\rangle = \left\langle {\Omega }_{II}\right|\mathop{\prod }\limits_{{\bf{x}},k}{w}_{II}({\bf{x}},k)\mathop{\prod }\limits_{{\bf{x}},k}{{\mathcal{U}}}_{II}^{{\mathcal{G}}}({\bf{x}},k) \mathop{\prod }\limits_{{\bf{x}}}{A}_{II}({\bf{x}})\left|{\Omega }_{II}\right\rangle \left|{\Omega }_{{\rm{p}}}\right\rangle \end{array}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left|{\Omega }_{J}\right\rangle$$\end{document} is the Fock vacuum for the virtual modes used in each construction; \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left|{\Omega }_{{\rm{v}}}\right\rangle =\left|{\Omega }_{I}\right\rangle \left|{\Omega }_{II}\right\rangle$$\end{document} . Similarly, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left|{\Omega }_{{\rm{p}}}\right\rangle$$\end{document} is the Fock vacuum of the physical matter. The first two products are over all links while the last is over all sites. If we ignore the portion of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${U}_{J}^{{\mathcal{G}}}$$\end{document} which consists of an operator on the gauge fields, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\psi }_{I}({\mathcal{G}})$$\end{document} is a complex number while \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left|{\psi }_{II}({\mathcal{G}})\right\rangle$$\end{document} is a state of the matter fields.

Substituting into Equation (21), the full ansatz can be written

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{array}{rcl}\left|\psi \right\rangle = \mathop{\sum }\limits_{{\mathcal{G}}}\,\left\langle {\Omega }_{{\rm{v}}}\right|\mathop{\prod }\limits_{{\bf{x}},k}\omega ({\bf{x}},k)\mathop{\prod }\limits_{{\bf{x}},k}{{\mathcal{U}}}^{{\mathcal{G}}}({\bf{x}},k) \mathop{\prod }\limits_{{\bf{x}}}A({\bf{x}})\left|{\Omega }_{{\rm{v}}}\right\rangle \left|{\Omega }_{{\rm{p}}}\right\rangle \big|{\mathcal{G}}\big\rangle \end{array}$$\end{document}

which satisfies all that we wished from a state: guaranteed gauge invariance, an entanglement area law, and efficient representation for calculation of observables. Note that the only variational parameters of the state are contained in A by way of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal{T}}$$\end{document} .

We finish this section with a brief review of some properties of this state. We first note that for any particular gauge field configuration \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal{G}}$$\end{document} , the state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\psi }_{I}({\mathcal{G}})\left|{\psi }_{II}({\mathcal{G}})\right\rangle$$\end{document} is Gaussian. As a result, this fermionic state is fully characterized by its covariance matrix, which can be computed efficiently, and from which (together with the gauge field configuration) all observables can be calculated. We therefore satisfy the first requirement—computationally efficient calculation of observables.

It is known that local, gapped Hamiltonians in one dimension obey an entanglement area law—that is, the entanglement between subsystems grows proportionally to the size of the boundary between them, which in one dimension is simply a constant^30,31^. The analogous statement is conjectured rather than proven in two dimensions, but our ansatz satisfies this property automatically—the entanglement between subsystems can arise only due to the entanglement introduced by the projection operators w(x, k), and along any curve which divides the lattice into two, the number of such operators is proportional to the length of the curve^61^. This is illustrated in Fig. 3.Fig. 3A graphical representation of the state which shows that an entanglement area law is automatically obeyed.The physical matter on the lattice sites is shown in yellow, the virtual modes are shown in turquoise, and the gauge fields are not shown. The entanglement between the region highlighted in purple and the rest of the lattice is proportional to the number of links through which the boundary cuts, since the only source of entanglement is the projection onto maximally entangled states of the virtual modes on the same link (enclosed by the turquoise dashed lines).

Finally, our state is gauge invariant — it is an eigenstate of the operator defined in Equation (8). In particular, it satisfies Equation (9) (or equivalently Equation (17), which defines gauge-invariance for the transformed fermionic operators), though the ansatz is general enough to handle other sectors as well^61^.

Algorithm

We wish to study the system described in the Physical System subsection of the Methods section by finding its ground state. To do so, we employ a variational Monte Carlo procedure. The key insight is that the expectation value of any observable can be written as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\langle {\mathcal{O}}\right\rangle =\mathop{\sum }\limits_{{\mathcal{G}}}{F}_{{\mathcal{O}}}\left({\mathcal{G}}\right)p\left({\mathcal{G}}\right),$$\end{document}

where the summation is over all possible gauge field configurations (in the results presented here, we do not make use of gauge fixing). \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${F}_{{\mathcal{O}}}\left({\mathcal{G}}\right)$$\end{document} can be expressed in terms of the covariance matrices of the state, and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p({\mathcal{G}})=\frac{{\left|{\psi }_{I}({\mathcal{G}})\right|}^{2}\langle {\psi }_{II}({\mathcal{G}})| {\psi }_{II}({\mathcal{G}})\rangle }{{\sum }_{{\mathcal{G}}{\prime} }{\left|{\psi }_{I}({\mathcal{G}}{\prime} )\right|}^{2}\langle {\psi }_{II}({\mathcal{G}}{\prime} )| {\psi }_{II}({\mathcal{G}}{\prime} )\rangle }$$\end{document}

defines a valid probability distribution that can be used for Monte Carlo. The expressions for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${F}_{{\mathcal{O}}}\left({\mathcal{G}}\right)$$\end{document} for the electric and magnetic terms of the energy are given in ref. ^63^. The expressions for the interaction and mass terms are given in Supplementary Information B: Observables.

To find the ground state, the relevant observables are the gradients of the energy — in particular, the gradients with respect to the parameters that define a particular state, which are all contained in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal{T}}({\bf{x}})$$\end{document} in A(x). We start by initializing the state with random parameters (subject to the constraints of Supplementary Information A: Ansatz Details), evaluate the gradients with Monte Carlo, and then update the parameters using the BFGS optimization algorithm until convergence is reached.

It is also worth noting that in the Monte Carlo procedure, we must calculate observables for many gauge field configurations, but that this does not require calculating the covariance matrix from scratch. Instead local updates can be applied at intermediate steps in the calculation depending on which gauge fields were modified^61^. This further improves the computational efficiency of this approach.

Results

In the previous sections we described our ansatz state and minimization procedure for finding the ground state of the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{2}$$\end{document} Kogut-Susskind Hamiltonian shown in Equation (7). Some details regarding the configuration of Monte Carlo and the minimization are given in Supplementary Information D: Computation Details. As a first check of our method, we compare our results with results found by exactly diagonalizing the Hamiltonian. This could only be done (under practical time and memory constraints) for a 2 × 2 and 4 × 4 lattice. Note that due to the staggering of matter, the extent of the lattice in each dimension must be an even number. Figure 4 shows a scan of the electric (and magnetic) couplings, with the corresponding ground state energy as well as each term of the energy—compared with exact diagonalization results.Fig. 4. The various energy observables for the 2 × 2 ground states.The dots show results given byour ansatz, showing their close match to the exact diagonalization results (solid lines) for each term in the Hamiltonian. The interaction coupling gI was fixed at gI = 1.0, and the mass coupling gM was fixed at gM = 1.0. The horizontal axis is the electric coupling gE = λ; the magnetic coupling is defined as gB = 1/λ. The results shown here do not rely on Monte Carlo, but rather on exact contraction --- for a 2 × 2 system, it is possible to simply iterate over all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${2}^{2{L}^{2}}=256$$\end{document} gauge configurations of the entire lattice. No error bars are shown as no sampling error is present and we have not quantified numerical error due to the minimization procedure.

The most important result from Fig. 4 is that our ansatz not only finds the ground state energy, but also correctly describes the physics of the ground state, as seen from the close agreement of each term in the energy with the exact results. Note that since the state is found using a variational procedure it provides an upper bound on the total energy of the true ground state, but there is no such guarantee on each term.

The couplings in the Hamiltonian (Equation (7)) define a three dimensional parameter space (since an overall energy scaling is irrelevant, we can fix one of the couplings). A two dimensional slice of this parameter space is shown in Fig. 5, with the results again calculated using exact contraction for a 2 × 2 system. The figure shows the value of 1 × 1 Wilson loops. Though Wilson loops are not an order parameter when fermionic matter is included, the transition between high and low values perhaps signals a phase transition between confined and unconfined phases; we leave an investigation of phases for future work.Fig. 5A parameter diagram for a 2 × 2 system with massless fermions (gM = 0) coupled to gauge fields.The horizontal axis is the electric coupling gE = λ (the magnetic coupling is given by 1/λ), and the vertical axis is the interaction coupling gI. The colors indicate the value of 1 × 1 Wilson loops for the ground state of the Hamiltonian with the specified couplings.

Figure 6 shows the ground state energy for lattice sizes 2 × 2, 4 × 4, and 6 × 6 as a function of all of the couplings. The results for the 2 × 2 were calculated using exact contraction, while the 4 × 4 and 6 × 6 results used Monte Carlo. The top row shows massless fermions (gM = 0) while the bottom row shows massive fermions (gM = 1.0). The results in the upper left pane (2 × 2 lattice with massless fermions) correspond to horizontal slices of Fig. 5; the results in the lower left pane (2 × 2 lattice with massive fermions) with gI = 1.0 are the same as the total energy shown in Fig. 4. The 4 × 4 exact-diagonalization data (which required  > 100 GB of memory for each point) makes it clear that the ansatz has some difficulty with the optimization at the peak energy, though with greater computation time, the convergence could likely be improved.Fig. 6. The ground state energies found for various lattice sizes and couplings.Moving from left to right, panels (a-b) show a 2 × 2 lattice (computed via exact contraction), panels (c-d) show a 4 × 4 lattice (computed with Monte Carlo), and panels (e-f) show a 6 × 6 lattice (computed with Monte Carlo). Where available, the results computed via exact diagonalization are shown as solid lines. The top row (panels a, c, e) is systems with gM = 0 (massless fermions), while the bottom row (panels b, d, f) gives results for gM = 1. The interaction coupling gI is shown in the legend, while λ on the x-axis defines the electric and magnetic couplings (gE = λ and gB = 1/λ). The error bars shown in the legend (which are generally too small to be seen on the graph) are the errors due to Monte Carlo sampling; error due to imperfect minimization is not shown.

Finally, we demonstrate that our method can be extended to larger lattice sizes, which are beyond the reach of exact diagonalization, and that we capture the physics of Wilson loops for such systems. Figure 7 shows Wilson loops for 4 × 4 and 6 × 6 systems. Here too we see a sharp change in the values of Wilson loops across the transition in ground state energies shown in Fig. 6.Fig. 7. The expectation values of 2 × 2 Wilson loops for 4 × 4 and 6 × 6 systems.Though Wilson loops are not an order parameter for theories including matter, the sharp transition is noteworthy. The interaction coupling gI is shown in the legend, while λ on the x-axis defines the electric and magnetic couplings (gE = λ and gB = 1/λ). The error bars shown in the legend (which are generally too small to be seen on the graph) are the errors due to Monte Carlo sampling.

Discussion

This paper demonstrated the use of GGPEPS for studying a lattice gauge theory with dynamical fermions. Previous work has shown that GGPEPS can capture the low energy physics of pure-gauge \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{2}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{3}$$\end{document} theories. The key advance demonstrated here is the inclusion of fermionic matter — in addition to gauge fields — in the model. This is a significant further step on the path to simulating lattice gauge theories in more than one spatial dimension, in cases where they cannot be efficiently simulated using standard Monte Carlo approaches due to the sign problem. Our results agree with exact results where they can be compared, and show reasonable behavior for larger system sizes where they cannot.

Our ansatz is a superposition of Gaussian states coupled to gauge configurations, which allows us to use the covariance matrix formalism in order to efficiently calculate observables. Since it is a superposition involving Gaussian states, our final ansatz \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left|\psi \right\rangle$$\end{document} is not itself Gaussian. This has proven sufficient for generating the numerical results demonstrated here, including the interaction between gauge fields and matter; more complex LGTs may require extending the state search over superpositions of these \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left|\psi \right\rangle$$\end{document} , each with different parameters in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathcal{T}}$$\end{document} , as discussed in^70^. Together with Monte Carlo integration over the gauge configurations, our approach allows for a variational ground state search which demonstrates the ability of the ansatz to capture the low-energy physics of a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\mathbb{Z}}}_{2}$$\end{document} lattice gauge theory with matter. In total, this method achieves our three aims: efficient representation of and computation with the state, guaranteed gauge invariance, and the entanglement area law expected for the system under consideration.

The next step in applying GGPEPS to LGTs is their application to a model that suffers from the sign problem, such as one with several flavors of fermions with different chemical potentials. Recent updates to our simulation code have significantly improved its efficiency, which will allow for increasing the size and complexity of the systems that can be studied using this approach, as well as greater investigation of the numerical performance of the algorithms as a function of ansatz settings. This will pave the way for studying theories with other gauge groups, culminating with the study of quantum chromodynamics in 3+1d.

Supplementary information

Supplementary Material