A variational Monte Carlo algorithm for lattice gauge theories with continuous gauge groups: a study of (2+1)-dimensional compact QED with dynamical fermions at finite density

Lattice gauge theories coupled to fermionic matter account for many interesting phenomena in both high energy physics and condensed matter physics. Certain regimes, e.g. at finite fermion density, are difficult to simulate with traditional Monte Carlo algorithms due to the so-called sign-problem. We present a variational, sign-problem-free Monte Carlo method for lattice gauge theories with continuous gauge groups and apply it to (2+1)-dimensional compact QED with dynamical fermions at finite density. The variational ansatz is formulated in the full gauge field basis, i.e. without having to resort to truncation schemes for the $U(1)$ gauge field Hilbert space. The ansatz consists of two parts: first, a pure gauge part based on Jastrow-type ansatz states (which can be connected to certain neural-network ansatz states) and secondly, on a fermionic part based on gauge-field dependent fermionic Gaussian states. These are designed in such a way that the gauge field integral over all fermionic Gaussian states is gauge-invariant and at the same time still efficiently tractable. To ensure the validity of the method we benchmark the pure gauge part of the ansatz against another variational method and the full ansatz against an existing Monte Carlo simulation where the sign-problem is absent. Moreover, in limiting cases where the exact ground state is known we show that our ansatz is able to capture this behavior. Finally, we study a sign-problem affected regime by probing density-induced phase transitions.


I. INTRODUCTION
Gauge theories play a prominent role in different areas of physics.In high-energy physics, the standard model of particle physics, a gauge theory, describes three of the four fundamental forces in nature.At high energy scales its interactions can be treated perturbatively, however, at lower energies this approach fails and non-perturbative techniques are required [1,2].This naturally gives rise to lattice gauge theories as they are non-perturbative, gauge-invariant regularizations of quantum field theories [3,4].In condensed matter, lattice gauge theories emerge as low-energy effective theories of strongly correlated electron systems, e.g.quantum spin liquids or high-temperature superconductors [5,6].
Much progress has been made in studying lattice gauge theories, both from the high-energy physics as well as from the condensed matter side, in particular using Euclidean Monte Carlo simulations [7].Nevertheless, certain regimes are difficult to access within this framework as fermionic theories at finite density or with an odd number of fermion flavors may suffer from the sign problem [8] and real-time dynamics are difficult to compute as Monte Carlo algorithms are usually formulated in Euclidean spacetime.
In recent years, several approaches to this problem have received attention: a prominent example is quantum simulation where it was shown that lattice gauge theory Hamiltonians can be realized in quantum devices (e.g.ultracold atoms, trapped ions or superconducting qubits) [9,10].The implementation of quantum simulators has been demonstrated in one dimension using trapped ions and ultracold atoms [11][12][13][14][15].In two and more spatial dimensions the situation becomes more challenging, in particular due to appearance of magnetic interactions, leading to four-body plaquette terms on the lattice.There have been proposals on how to overcome this problem in quantum simulators (either by employing a digital [16][17][18][19][20][21] or an analog simulation scheme [22,23]) but so far they have not been realized in experiments.
Another significant approach is based on variational ansatz states which can capture the relevant physics of the theory but at the same time can be evaluated efficiently.For lattice gauge theories these states either have to respect the local gauge symmetries or one has to find a reformulation of the theory in terms of gaugeinvariant variables (at the cost of more complicated interactions) [24][25][26] such that there is a larger freedom in choosing variational states.One class of ansatz states are tensor networks whose one-dimensional version, matrix product states (MPS), have been successfully applied to (1+1)-dimensional Abelian and non-Abelian lattice gauge theories [27][28][29][30][31][32][33][34][35][36][37], enabling the study of finite chemical potential scenarios and out-of-equilibrium dynamics which are not accessible in Monte Carlo simulations of Euclidean lattice gauge theory.In higher dimensions, tensor network methods have been applied to lattice gauge theories with a finite-dimensional gauge field Hilbert space (either by working with quantum link formulations [38][39][40][41] or using a discrete gauge group [42][43][44]).Other types of ansatz states can be formulated in the infinite-dimensional Hilbert space of continuous gauge groups.These include periodic Gaussian states [45] (generalizations of Gaussian states that take into account the compactness of the gauge group) or neural network-based ansatz states [46,47].
A particularly interesting model in higher dimensions is (2 + 1)-dimensional compact QED (cQED 3 ) with dynamical charges, both in the context of high-energy physics and condensed matter physics.From the highenergy perspective it is interesting since it is the simplest theory to discuss confinement and chiral symmetry breaking, also quintessential to our understanding of quantum chromodyanmics (QCD).Already without dynamical fermions, the theory has non-trivial interactions due to the appearance of four-body magnetic terms.It is known to confine for all couplings [48].However, upon the inclusion of dynamical fermions, the situation is less clear since the dynamical fermionic matter might lead to deconfinement.These phenomena are also of high relevance in condensed matter since many low-energy effective theories of two-dimensional strongly coupled electron systems can be described by massless dynamical fermions coupled to a compact U (1) gauge field.Compared to Z 2 lattice gauge theories where it was shown that sign-problem-free Monte Carlo simulations could be performed for an even number of fermion flavors even at non-zero density [49], for the U (1)-theory the sign-problem is only absent for an even number of fermion flavors at half-filling [50].Based on the above considerations, in this work, we introduce a variational method that can access the sign-problem affected regimes of cQED 3 with dynamical charges without truncating the U (1) gauge field Hilbert space.The ansatz is based on a combination of a pure gauge part containing the selfinteractions of the gauge field and a fermionic part which describes the dynamics of the matter degrees of freedom with the gauge field.The pure gauge part is a Jastrow wave function constructed out of gauge-invariant plaquette variables (its form can be connected to certain neural network quantum states [51]).The choice of ansatz is motivated by an earlier proposal [45] which could approximate ground states and real-time dynamics in cQED 3 with static charges.The fermionic part is an infinite superposition of gauge-field dependent fermionic Gaussian states which are parametrized in such a way that the resulting state is gauge-invariant.Note that the parametrization is done in a way that the number of parameters only scales polynomially with system size.In a similar fashion to neural-network quantum states, expectation values are obtained using Monte Carlo sampling.The optimization of variational states is done via stochastic reconfiguration [52].
In order to verify the capabilities of the variational method we first demonstrate also numerically that gauge invariance is indeed preserved.In a second step it is shown that the ansatz is exact in all limiting cases.This includes the weak-coupling limit (g 2 → 0) where the ground state is known to be a π-flux state [53] and the strong-coupling limit (g 2 → ∞) where the electric energy dominates and one obtains an effective fermionic theory.Moreover, we benchmark our ansatz against other meth-ods: first, we compare for the pure gauge theory, i.e. compact QED without fermions, the ground state energy with another variational method [45] and see agreement over the whole coupling region.These results were recently confirmed in another variational study based on neural-network states [46].For compact QED with dynamical fermions, we compare with a Euclidean Monte Carlo study at zero chemical potential and two fermion flavors where the sign-problem is absent [50].We compute the flux energy per plaquette which agrees with ref. [50] over the whole coupling region.We also compare fermionic correlations quantifying the degree of antiferromagnetic order in the ground state.By extrapolating this quantity to the thermodynamic limit it is shown that antiferromagnetic order only persists down to a coupling of g 2 c,∞ = 0.15 (2) which is in qualitative agreement with ref. [50] although our extrapolated value for the transition is lower.To demonstrate our method in a signproblem affected regime we study density-induced phase transition for two fermion flavors at non-zero chemical potential, similar to a tensor network study in one dimension [33].We consider both the case of massless and massive staggered fermions and see qualitatively similar phenomena as in ref. [33].
The manuscript is structured as follows.In section II, we introduce the model, cQED 3 with dynamical fermions.In section III, we describe the variational Monte Carlo method including the gauge-invariant construction of the variational state, the numerical evaluation with Monte Carlo techniques and the adaptation of variational parameters.In section IV, our ansatz is benchmarked against limiting cases of the model where the ground state is known and other numerical methods.
In section V, we study a sign-problem affected regime by investigating density induced phase transition for two fermion flavors at non-zero chemical potential.In section VI, we summarize and conclude.

II. THE MODEL: (2 + 1)-DIMENSIONAL
COMPACT QED WITH DYNAMICAL FERMIONS We study (2 + 1)-dimensional compact quantum electrodynamics (cQED 3 ) coupled to dynamical fermions.The model is defined on an L × L square lattice with periodic boundary conditions.We work with staggered fermions [54] which are suitable for studying chiral symmetry breaking.The fermions can appear in several species α which can be subject to different chemical potentials µ α (in some scenarios they are also given a mass FIG. 1. Naming conventions on the periodic lattice: the gauge field degrees of freedom θx,i, Ex,i (blue) reside on the links x, i while the fermionic degrees of freedom ψ † x,α (red), which can come in several species α, are located on the sites x.The circular arrows on the plaquettes denote the plaquette variables θp.The global loops θi wind around the axis given by ei and are illustrated by blue lines.

m). The Hamiltonian reads
where ψ x,α denotes the fermionic annihilation operator for site x and species α.The gauge field operator θx,i and the electric field operator Êx,i fulfill the canonical commutation relations, [ θx,i , Êy,j ] = iδ ij δ x,y .Accordingly, the gauge field on a link can be represented either by an integer-valued electric field variable, Êx,i or by an element of the U (1) gauge group, θx,i |θ x,i = θ x,i |θ x,i (θ x,i ∈ [0, 2π)).We will mostly use the group element representation throughout the manuscript.In this representation, the electric field operator has the form Êx,i = −i∂/∂θ x,i .The plaquette operator θp = θx,1 + θx+e1,2 − θx+e2,1 − θx,2 is the clockwise summation of link operators around plaquette p where x is the site at the bottom left corner.The labelling conventions are illustrated in Fig. 1.The magnetic coupling g mag is usually chosen to be where the staggered charge operator Qstag,x is defined as Physical states |phys must be eigenstates of all Gauss law operators Ĝx |phys = q x |phys ∀x (4) where eigenvalues q x correspond to different static charge configurations.

III. THE VARIATIONAL METHOD
Since our method is based on variational Monte Carlo, we will explain it in several steps (sketched in Fig. 2): we first discuss the state construction and motivate the choice of our ansatz.In the second step, we explain the evaluation procedure of our ansatz based on Monte Carlo sampling.In the third and final step, we discuss the adaptation of variational parameters based on stochastic reconfiguration.

A. State construction
We construct our ansatz state in the gauge field basis where states are characterized by all U (1) gauge link variables θ x,i , |{θ x,i } ≡ ⊗ x,i |θ x,i .A general gauge field state can then be defined by where Ψ G ({θ x,i }) is a function over all gauge link variables θ x,i .To extend the above to an arbitrary state |Ψ of the cQED 3 model introduced in section II (i.e.including fermions) we need to specify a fermionic Fock state |Ψ F ({θ x,i }) for every gauge field configuration {θ x,i }: where we abbreviated for ease of notation the gauge field configurations {θ x,i } as θ and the measure as Dθ = x,i 2π 0 dθ x,i .This notation will be used throughout the article.
One should note that an arbitrary state |Ψ as given above is a priori not gauge-invariant.Thus, the gauge invariance condition for physical states in eq. ( 4) severely restricts the possible choices for Ψ G (θ) and |Ψ F (θ) .
The state |Ψ defined above is completely general.From now on we will use the form of |Ψ as the basis to construct our variational ansatz state which will be defined by specifying the pure gauge part Ψ G (θ) and the fermionic ansatz |Ψ F (θ) .
Intuitively, the role of Ψ G (θ) and |Ψ F (θ) in our construction can be motivated as follows: Ψ G (θ) is designed to approximate the ground state of the pure gauge model H KS ≡ H B + H E (the Kogut-Susskind Hamiltonian [4]) whereas |Ψ F (θ) is designed to approximate the low-energy physics of the fermionic Hamiltonian H fer ≡ H E + H GM + H M which neglects the self-interactions of the gauge field.

The pure gauge part of the ansatz
In this section we motivate and describe the pure gauge part Ψ G (θ) of our variational state.In earlier work [45] it was shown that is a good ansatz for the ground state of compact QED with static charges.It is a Gaussian in the plaquette variables θ p that is made periodic by an infinite sum over the integer-valued variables N p (α pp are variational parameters).The periodicity is important to account for the compactness of the U (1) gauge field.
Here, we would like to find an ansatz which has a similar expressive power as the states above but at the same time is suitable for a variational Monte Carlo simulation directly in θ (without resorting to the sums above).A useful hint is given by the Villain approximation [55] which states for γ → ∞.Therefore, a suitable ansatz state could be e − pp cos(θp)α p,p cos(θ p )+ p βp cos(θp) (9) with the matrix α and the vector β being variational parameters.We will choose α and β to be real since we are interested here in low-energy properties.For the study of real-time dynamics (which will be investigated in a future work) we would choose the variational parameters to be complex.Apart from the cosine terms we can add sine terms, combine them in a vector b(θ) = (cos(θ p1 ), .., cos(θ p N ), sin(θ p1 ), .., sin(θ p N )) and generalize the above state to which will be the variational ansatz for the pure gauge field dynamics entering the full ansatz as in eq. ( 6).
In the case of periodic boundary conditions there are two inequivalent global non-contractible loops (inequivalent in the sense that they can not be transformed into each other by plaquette operations).We choose them to be θ 1 (winding around the lattice along the x 1 -axis), and θ 2 (along the x 2 -axis), respectively (see Fig. 1).We incorporate them in our ansatz by expanding the vector b(θ) by the entries cos(θ 1 ), cos(θ 2 ), sin(θ 1 ) and sin(θ 2 ).This is necessary because upon the coupling of compact QED to dynamical fermions, θ 1 and θ 2 become dynamical variables due to the appearance of gauge-matter interactions where the phase e iθx,i appears.If expressed in terms of gauge-invariant variables, it contains contributions from both plaquette variables θ p and the global loops θ 1 and θ 2 [25].For static charges, the magnetic Hamiltonian is the only term depending on the gauge field variables θ x,i which can be expressed entirely in terms of plaquette variables θ p such that the global loop variables only set different topological sectors (similar to the toric code).
For all our purposes it turned out that all variational parameters in α corresponding to the global loop variables were not relevant and that it was sufficient to only keep the global loop parameters in β variational.After imposing translational invariance we thus remained with 2N + 4 variational parameters for α and 6 variational parameters for β (with N = L 2 the number of lattice sites).
Since b(θ) contains only closed loops the gauge field part Ψ G (θ) as a function of b(θ) automatically preserves gauge invariance.

The fermionic part of the ansatz: gauge-invariant fermionic Gaussian states
The fermionic part of our variational ansatz |Ψ F (θ) is a generalization of fermionic Gaussian states that can incorporate interactions between fermions and gauge fields while preserving gauge invariance.The overall state Dθ |Ψ F (θ) |θ is an integral over all gauge field configurations where for every gauge field configuration θ we define a fermionic Gaussian state |Ψ F (θ) .The motivation for this construction is that the resulting state, an infinite superposition of Gaussian states, is a powerful ansatz state as it is clearly not Gaussian anymore and can capture correlations beyond the Gaussian realm.At the same time, we retain for every |Ψ F (θ) the properties of fermionic Gaussian states which allows us to compute part of the expectation values analytically.The number of variational parameters is shown to scale only polynomially in the system size and not exponentially as the number of gauge field configurations.
Recalling that every pure fermionic Gaussian state can be represented by a unitary operator U GS acting on some reference state |Ψ 0 [56], we carry out an analogous procedure for every gauge field configuration θ to construct gauge-field dependent fermionic Gaussian states as In our method, the reference state |Ψ 0 will be chosen as the ground state of H fer = H E + H GM + H M in the strong-coupling limit (g 2 → ∞), i.e. the regime where the electric term dominates so that electric field excitations are strongly suppressed.In the following we will refer to strong-and weak-coupling always w.r.t. the relative strength of the electric Hamiltonian H E (quantified by its coupling constant g 2 ).If we only consider H E in the strong-coupling limit, |Ψ 0 will be a Fock state where all fermions are fixed to certain sites (the exact form of the state will depend on the number of fermion flavors and the configuration of background charges).However, |Ψ 0 does not need to be Gaussian but one can also include perturbations induced by gauge-matter interactions H GM , e.g. it is known that for two fermion flavors at half-filling the strong-coupling ground state in second-order perturbation theory is the ground state of the Heisenberg model.It can be shown that its properties can be incorporated in the reference state |Ψ 0 which can even be kept variational (for details see Appendix C).
For simplicity of the discussion, we will assume in the following a Gaussian reference state and only one flavor of staggered fermions (how the ansatz can be readily extended to multiple flavors is described in Appendix B).In the sector without background charges the reference state is chosen to be the Dirac state i.e. with all odd sites O occupied.
The Gaussian operator U GS (θ) acting on |Ψ 0 is defined as where ξ(θ) is dependent on the gauge-field and on the variational parameters.The gauge-field dependence has to be chosen in a way that respects gauge invariance.This can be achieved by defining ξ(θ) via the eigendecomposition of the gauge-matter Hamiltonian which can be written as with ψ x ≡ (ψ x1 , ..., ψ x N ) T a vector of all fermionic annihilation operators.The matrix h GM (θ) is hermitian and can be diagonalized for a specific gauge field config- with ξ containing the variational parameters.Note that ξ does not depend on the gauge field configuration and thus the number variational parameters scales quadratically with the system size (linearly for our choice of parametrization, see Appendix B).Putting everything together, the fermionic part of the ansatz for one fermion flavor takes the form and the whole variational ansatz state |Ψ is thus fully defined according to eq. ( 6).Gauge invariance of |Ψ follows from the fact that H GM and its eigenstates are gauge-invariant since the construction of |Ψ F (θ) given in eq. ( 16) is formulated in terms of these eigenstates.
Since the gauge invariance condition in eq. ( 4) is local in θ, every realization of the state in a Monte Carlo simulation will be gauge-invariant, i.e. even with an imperfect sampling algorithm the unphysical part of the Hilbert space is never accessed.
The motivation for the choice of ansatz above is on the one hand that it ensures gauge invariance but more importantly, by choosing the matrix ξ appropriately, the occupation of the eigenstates of H GM can be tuned which allows to obtain good ground state approximations even in regimes where strong gauge field fluctuations are present.This has to be seen in contrast to mean-field descriptions where a certain gauge field pattern is fixed and the resulting fermionic theory is studied.The latter has the problem, which is particularly relevant in the study of quantum spin liquid states (where the lattice gauge theory emerges as an effective low-energy description), that it often remains unclear whether the spin liquid state is stable against gauge field fluctuations [5].
The cost of working with the ansatz is that the eigendecomposition of h GM (θ) needs to be carried out at every measurement step of the Monte Carlo algorithm.However, h GM (θ) is a hermitian N ×N matrix where N is the number of lattice sites such that the cost is O(N 3 ) which can be done efficiently.Note that the number of fermion flavors does not enter as the gauge-matter interaction is the same for all flavors.
The fermionic ansatz state in eq. ( 16) is normalized since the Gaussian operator acting on the Dirac vacuum is unitary.This is beneficial for the variational Monte Carlo simulation since it will not contribute to the probability distribution that needs to be sampled.Thus, no sampling problems related to fermion determinants can occur in this method as opposed to action-based Monte Carlo algorithms.
So far we have not specified the matrix ξ in eq. ( 16) containing the fermionic variational parameters.For that we consider the eigendecomposition of h GM (θ), denoted as h GM (θ) |w i (θ) = λ i (θ) |w i (θ) , i ∈ {1, .., N }.Assuming an L × L lattice with L even, the spectrum of h GM (θ) is symmetric around zero, i.e. we have N/2 pairs of eigenvectors |w k+ (θ) and |w k− (θ) (k ∈ {1, .., N/2}) such that |w k+ (θ) corresponds to the eigenvalue +λ k (θ) and |w k− (θ) to the eigenvalue −λ k (θ).A useful feature of these pairs is their structure in the position basis as they can be written as two vectors |w ke (θ) , |w ko (θ) which are residing exclusively on even (respectively odd) lattice sites: This allows us to write the strong-coupling state in eq. ( 12), where the fermions occupy all odd sites, as a product over all pairs k where in each pair the odd superpositon is occupied, ).The purpose of ξ in eq. ( 16) is then to smoothly transform this equal superposition of |w k+ (θ) and |w k− (θ) into a state where all |w k− (θ) are occupied, corresponding to the ground state of H GM .Thus, ξ allows us to transform smoothly from the strong-coupling ground state to the weak-coupling ground state.For more details on ξ and the specific choice of parametrization see Appendix B.

B. Evaluating expectation values
In this section we describe how Monte Carlo techniques can be used to compute various expectation values for the variational ansatz presented in the previous section.Throughout the following discussion the variational parameters are kept fixed, their adaptation will be discussed in the next section.
For the computation of an observable O with the full ansatz |Ψ from eq. ( 6) we obtain where |Ψ F (θ) is absent in the norm since it is already normalized by construction (see eq. ( 16)) so that the probability distribution p(θ) depends only on Ψ G (θ).The fermionic part of the ansatz thus only appears in the numerator for the evaluation of O which is carried out analytically and only the remaining expression O loc (θ) is sampled in a Monte Carlo simulation.
We split the calculation of O loc (θ) in two parts: since O is a priori not diagonal in θ (e.g.all electric observables involve derivatives w.r.t.θ) we first compute the action of O on our ansatz |Ψ which gives rise to an expression O fer (θ) that is diagonal in θ but might still contain fermionic operators (e.g.due to derivatives of the fermionic ansatz |Ψ F (θ) ): which is now a real-valued function that can be readily sampled in a Monte Carlo simulation.The probability distribution p(θ) according to which we need to sample is only dependent on the gauge part Ψ G (θ) defined in eq.(10): The method described above has to be contrasted with usual variational Monte Carlo methods [57] where the whole trial wavefunction contributes to the probability distribution and the local quantities O loc (θ) do not involve taking expectation values w.r.t.some part of the ansatz.
Having discussed the general procedure, the computation of observables can be divided into three groups by level of difficulty: the first group consists of observables that are not diagonal in θ (all electric quantities such as H E ) and thus first need to be brought into a diagonal form O fer (θ).These observables are the most involved.
The second group of observables are already of that form but since O fer (θ) is still a fermionic operator it needs to be evaluated w.r.t.|Ψ F (θ) to obtain O loc (θ) (e.g.H GM or H M ).The third group of observables is already of the form O loc (θ) and can be readily sampled in a Monte Carlo simulation (e.g.H B ).
Two things need to be shown to demonstrate that our variational ansatz can be used efficiently: first, the efficient computation of O loc (θ) and secondly, efficient sampling of the probability distribution p(θ).Thus, in the following, we first show exemplary for the Hamiltonian how O loc (θ) is derived, i.e. we compute the local energy H loc (θ).In a second step, we explain the Monte Carlo simulation, in particular how samples from p(θ) are generated using Metropolis algorithm.

Computation of the local energy H loc (θ)
The electric Hamiltonian H E is the only term of the Hamiltonian defined in eq. ( 1) that is not diagonal in θ (the most difficult type of observable to compute, as discussed above).We thus focus on H E and discuss other terms briefly at the end of this section.
The electric Hamiltonian corresponds to second order derivatives in the gauge field variables θ x,i .Since our ansatz consists of a fermionic part |Ψ F (θ) and a pure gauge part Ψ G (θ), the electric energy has a solely fermionic contribution, a pure gauge contribution and a crossterm between the two, denoted as: We start by considering H E gg , the part originating from taking twice the derivative of Ψ G (θ) whose construction is based on the vector b(θ) (see eq. ( 10)).Hence, we need to compute the derivative of b(θ) with respect to θ x,i which gives rise to the vector b x,i (θ) = δ p,(x,i) (− sin(θ p1 ), .., − sin(θ p N ), cos(θ p1 ), .., cos(θ p N )) with where (x, i) ∈ p clockwise (anti-clockwise) means that the link (x, i) is contained in the plaquette p and the orientation of the link is parallel (anti-parallel) to the orientation of the plaquette.For periodic boundary conditions we have the additional entries cos(θ j ) and sin(θ j ) in b(θ) corresponding to the global loops θ 1 and θ 2 .They give rise to the derivatives − sin(θ j ) and cos(θ j ) if (x, i) lies on the x j -axis and otherwise zero.The electric energy of the pure gauge part and the corresponding local quantity H E,gg,loc (θ) is then derived as Dθ e −S(θ) ≡ DθH E,gg,loc (θ) p(θ) with the probability distribution p(θ) and S(θ) = b T (θ)αb(θ) + 2β T b(θ), both defined in eq. ( 22).The part of the electric Hamiltonian acting only on Ψ G (θ) can therefore be written in a simple diagonal form in the gauge field basis.It is more difficult to compute H E,f f,loc (θ), i.e. the local quantity corresponding to derivatives of the fermionic ansatz |Ψ F (θ) .As discussed earlier, we first derive an expression H E,f f,fer (θ) that will be diagonal in θ but still contains fermionic operators (see Appendix A for details): The form of f x,i (θ) above is for a general gauge-field dependent fermionic Gaussian state characterized by some ξ(θ).To get an expression explicitly diagonal in θ we insert our ansatz ξ(θ) = V (θ) ξV † (θ) defined in eq. ( 15) which is based on the eigendecomposition of the gaugematter Hamiltonian, h GM (θ) = V (θ)Λ(θ)V † (θ).We obtain (see Appendix A for the derivation): We can find an explicit expression for α x,i (θ) which amounts to finding the derivatives of the eigenvectors of h GM (θ): where λ i (θ) are the eigenvalues of h GM (θ).The final expression for H E,f f,fer (θ) is thus diagonal in θ but still a quartic fermionic operator.This form of the electric Hamiltonian intuitively illustrates that the gauge field mediates interactions between the fermions.
In the following we want to evaluate these fermionic interactions w.r.t. the fermionic state |Ψ F (θ) as shown in the last row in eq. ( 27) to compute the local electric energy H E,f f,loc (θ) that can then be measured in our Monte Carlo simulation.
As a side note we want to mention that |Ψ F (θ) in its general form defined in eq. ( 11) does not need to be Gaussian as one can also choose a Non-Guassian reference state |Ψ 0 .This might be useful if one is particularly interested in the strong-coupling regime (from the high-energy physics perspective one is usually interested in the weak-coupling region where the continuum limit is located).In the strong-coupling regime the electric field is strongly suppressed and the Hilbert space effectively reduces to a fermionic Fock space.Such models can be tackled by other many-body methods (e.g.tensor networks) which are not suitable for lattice gauge theories with infinte-dimensional local Hilbert spaces.One could combine our ansatz with such methods by carrying out the unitary transformation given by the fermionic Gaussian operator U GS (θ) (acting on top of |Ψ 0 ) so that the remaining expression can be evaluated w.r.t. the reference state |Ψ 0 whose fermionic correlation functions could be computed with another method (in Appendix C we demonstrate this for two fermion flavors at half-filling where the effective model is the Heisenberg model).
If we focus, however, on the case of one fermion flavor and the Gaussian reference state |D as defined in eq. ( 16), we need to evaluate a fermionic Gaussian state for every gauge field configuration θ.The fermionic expectation values in eq. ( 27) can then be computed as where Γ ψψ † (θ) = V (θ) ΓV (θ) † is the covariance matrix of the Gaussian state |Ψ F (θ) and Γ = e i ξ V (θ) † Γ 0 V (θ)e −i ξ with Γ 0 the covariance matrix of the reference state |D and ξ containing the variational parameters (see eq. ( 15)).Inserting the expectation values above in eq. ( 27) gives H E,f f,loc (θ).
The last remaining part of the electric Hamiltonian, the crossterm H E f g , involves a quadratic expression in the fermions coming from |Ψ F (θ) and a derivative in b(θ) coming from Ψ G (θ) and is thus easier to compute than the quartic expressions in the pure fermionic contribution (for the explicit form see Appendix A).
Other parts of the Hamiltonian are easier to evaluate since they are already diagonal in the gauge field basis.For the sake of completeness we will provide them here briefly.First, the magnetic part which is directly suitable for Monte Carlo sampling: The gauge-matter interactions are already diagonal and only quadratic in the fermions: ≡ DθH GM,loc (θ)p(θ) (33) where the quadratic expressions in the fermions are evaluated in analogy to the electric part of the Hamiltonian.In the same fashion are other purely fermionic parts evaluated such as the mass term H M .
In terms of computational cost the local electric energy H E,loc (θ) is the most difficult part to evaluate.Naively, one expects the required number of operations for evaluating it to be O(N 4 ) (N the number of lattice sites) but with the chosen parametrization of ξ it can be shown to be O(N 3 ) (see Appendix B).

Monte Carlo algorithm
In the following we show how to efficiently evaluate an observable O with our ansatz |Ψ in a Monte Carlo simulation given an expression for O loc (θ).The expectation value of O is computed as an average over N samples θ i drawn from the probability distribution p(θ): The samples θ i are generated by a Markov chain θ 1 → . . .→ θ i → . . .→ θ N using Metropolis algorithm [58].One iteration in this procedure, i.e. θ i → θ i+1 , is described as follows: starting from θ i a new configuration θ is proposed according to some update scheme.In our case this involves sweeping through every link of the lattice and performing local updates on the gauge variables θ x,i .At the same time, we also perform global updates to switch between different monopole-like configurations which is hard to achieve with local updates (for details on the update scheme see Appendix D).Recalling from eq. ( 22) the form p(θ) ∼ e −S(θ) of our probability distribution, we compute the transition probability p(θ → θ ) = e −S(θ ) /e −S(θi) = e −∆S .In the acceptance step, a random number u between zero and one is generated and the new configuration is accepted if e −∆S ≥ u, i.e. θ i+1 = θ .Otherwise, the configuration θ is rejected and θ i+1 = θ i .In the first phase of the Monte Carlo simulation (the warm-up phase) these iterations are performed to equilibrate the system (i.e.reach configurations with sufficiently low weight S(θ)) and only after that the configurations θ i are used to compute the expectation value in eq.(34).
The numerical cost of performing Metropolis algorithm depends on computing the transition probability between the old configuration θ and the proposed configuration θ .For local updates they differ only in a single link variable θ x,i , respectively two plaquette variables θ p .The vector b(θ), constructed out of sin(θ p ) and cos(θ p ), is thus changed in four places.Since S(θ) is bilinear in b(θ), the cost of computing ∆S is of order O(N ) where N is the number of lattice sites.Sweeping through the lattice with this procedure is thus of order O(N 2 ).For the global updates the transition probability requires O(N 2 ) operations but is only performed O(1) times so that the cost of a full update is O(N 2 ).
Having such a low cost for updates has several advantages: we can perform multiple local and global updates to further decorrelate expensive measurements.The acceptance probability in our simulations stays on a high level throughout the whole coupling region (see Appendix D).Moreover, if we parallelize the Monte Carlo simulation with multiple runners there is practically no overhead due to the warm-up phase.

C. Adaption of variational parameters
In the last section we described the evaluation with our Monte Carlo scheme for a fixed set of variational parameters.To study ground states and dynamical phenomena we need to adjust the variational parameters accordingly.Here, we focus on the study of ground state properties but the discussion can be readily extended to time-evolution phenomena as we use an imaginary timeevoultion procedure (called stochastic reconfiguration in the variational Monte Carlo language [52]) to find the optimal set of parameters.We project the equations of motion onto the tangent plane of our variational manifold.For every variational parameter γ i , either fermionic (in ξ) or pure gauge (in α and β) we define a corresponding tangent vector |Ψ i ≡ P Ψ (∂ γi |Ψ ) where P Ψ ensures orthogonality to |Ψ : All tangent vectors in our ansatz are linearly independent which allows to invert the Gram matrix This can be intuitively explained by considering the different types of tangent vectors: the ones corresponding to the fermionic parameters are related to the singleparticle eigenstates of the gauge-matter Hamiltonian and are therefore orthogonal.The tangent vectors corresponding to the pure gauge part are quadratic (for α) or linear (for β) in the entries of the vector b(θ) which are related to the different plaquette variables θ p , thus leading to linearly independent tangent vectors.The imaginary time evolution of the variational parameters can then be expressed in the following way: with E ≡ Ψ|H|Ψ Ψ|Ψ the variational energy (whose evaluation was described in the previous section) and γ ≡ ∂γ ∂τ .The gradient of the variational energy and the Gram matrix need to be measured in a Monte Carlo simulation.The cost of both can be shown to scale in the same way as the cost of computing the variational energy (see Appendix E).Summarizing, the computational complexity of our variational Monte Carlo algorithm is O(N 2 ) for the update procedure and O(N 3 ) for the measurement procedure, thus allowing for an efficient implementation.

IV. BENCHMARKING OF THE VARIATIONAL METHOD
Now, we have all the ingredients to apply our variational method: We constructed a gauge-invariant state and showed how it can be efficiently evaluated for a fixed set of parameters using Monte Carlo sampling.Additionally, we have a scheme to adapt the parameters using stochastic reconfiguration.
In the following section, we investigate the validity of the variational method.It will be threefold: first, to confirm the analytical arguments about gauge invariance of the ansatz given in the previous section we will show numerically that our state is gauge-invariant up to machine precision.Secondly, we investigate different limiting cases of cQED 3 where the ground states are known.In the last part, we benchmark our results for the N f = 2 case at half-filling (in the sector of exactly one fermion per lattice site) with a recent Monte Carlo simulation [50].

A. Gauge invariance
To also show numerically that gauge invariance is manifest in our ansatz we compute the expectation value of the Gauss law operator G x (as defined in eq. ( 2)) for every site x and plot G x − q x for the whole lattice since this quantity needs to be zero for a physical, gauge-invariant state, see eq. ( 4).We choose different variational parameters, different lattice sizes and different sampling sizes but the violation of the Gauss law is always found to be of the order of machine precision, i.e.G x − q x 10 −16 .One such configuration for a very small sampling size of N = 10 and a system size L = 12 is illustrated in Fig. 3.

FIG. 3. Gauss law violation
Gx − qx for a 12 × 12 lattice at g 2 = 0.25, t = 1 and gmag = −1 with a random choice of variational parameters and a sampling size of N = 10.The Gauss law violation is of the order of machine precision even for a small sampling size, demonstrating that the ansatz is inherently gauge-invariant.

B. Limiting cases
It is useful to consider the limiting cases of compact QED with fermionic matter and convince ourselves that the ground state properties can be captured accurately by our method.In the following, we consider massless fermions without chemical potentials.
We first study the limit g 2 → 0 while keeping t and g mag fixed: it is well known that in this limit the gauge field forms a π-flux pattern and the fermions fill up the lower band at half-filling [53].A typical problem in meanfield theory is to investigate the stability of the π-flux pattern against gauge field fluctuations.This can be studied naturally in our ansatz by watching the parameter flow upon increasing the electric coupling constant g 2 .The π-flux state itself is naturally incorporated in our ansatz since we can fix the gauge field to a certain configuration by tuning the β-parameters to a very high value such that the constraint cos θ p = −1 is enforced for all plaquettes.In addition, since we have periodic boundary conditions, we also need to choose the optimal flux configuration for the global non-contractible loops which depends on the size of the lattice.To accomplish that it is important to have a global update in our update scheme since these global changes in the configuration can not be captured by only updating plaquettes locally.Finding the π-flux state is in general a useful test for our update scheme since the probability distribution needs to approximate a delta distribution for which a good update scheme is required.The fermionic part is obtained by tuning the variational parameters of the fermions in such a way that for all flux configurations the lower half of the band is oc- cupied (which corresponds to choosing all fermionic parameters ξ i = 1 as described in Appendix B).The result of our variational optimization is an accurate representation of the π-flux state with an average deviation on the order of 10 −8 from cos θ p = −1, respectively an average deviation on the order of 10 −4 from θ p = π (depicted in Fig. 4).
Next we consider the opposite limit to the π-flux state, the strong-coupling limit with large g 2 .In this limit, the electric energy dominates and some fluctuations are introduced in second-order perturbation theory by the gauge-matter Hamiltonian.For one fermion flavor this perturbation does not have a large effect and the ground state is described by a Gaussian state.For two fermion flavors, however, one can have correlated hopping processes which at half-filling give rise to the Heisenberg Hamiltonian.Both cases can be captured by construction in our ansatz since we design the fermionic part of the ansatz in such a way that our gauge-field dependent Gaussian operator acts on a strong-coupling reference state |Ψ 0 (see eq. ( 11)) and we can choose that reference state according to our needs.We can either choose a Gaussian state for |Ψ 0 or include more advanced methods, e.g. to approximate the Heisenberg ground state we can include spin wave theory in |Ψ 0 (see Appendix C).
We also benchmark for the limiting case that the gauge-matter interactions vanish (t = 0) so that fermions and gauge-field decouple and we obtain the standard pure gauge compact QED described by the Kogut-Susskind Hamiltonian H KS = H E + H B with g mag = 1 g 2 [4].We therefore only consider the pure gauge part of our ansatz (setting our fermionic variational parameters to Benchmark for pure gauge compact QED 3 of our ansatz (denoted by FG) against the variational method in ref. [45] based on complex periodic Gaussian states (denoted by CPG): we compare the ground state energy E0 for an 8 × 8 lattice over the whole coupling region and compute the relative error (see inset).
zero, ξ = 0).We benchmark our ansatz against a recent work [45] which has given good ground state and real-time dynamics of compact QED (the results were recently confirmed by another variational study [46]).We compare the ground state energy of both methods for an L = 8 × 8 for the whole coupling region of g 2 (since g 2 is the only coupling constant in pure gauge compact QED).We find that our results agree very well for the whole coupling region (with a maximal difference of half a percent) while our method performs a tiny bit better at large couplings where the method in ref. [45] gives minimally better results for small g 2 .The benchmark is illustrated in Fig. 5.

C. Benchmark against Euclidean Monte Carlo
Benchmarking for cQED 3 including dynamical fermions is in general difficult since in most scenarios a sign-problem occurs so that no Euclidean Monte Carlo studies exist.However, it was shown to be absent for an even fermion number at zero chemical potential [50].This was exploited in order to perform determinantal Monte Carlo simulations.Thus, it is natural to compare our ansatz with the Monte Carlo simulations for the case of N f = 2 fermionic species at half-filling.The analysis in ref. [50] revolves around the question of whether a confinement-deconfinement transition takes place and what the nature of this phase transitions is.
We fix the magnetic coupling and the gauge-matter coupling to g mag = −1 and t = 1 and will mostly vary the electric coupling g 2 .The first observable that is compared is the flux energy per plaquette cos(θ p ) averaged over the whole lattice.Our results are shown in Fig. 6.We see agreement over the whole coupling region of g 2 with Fig. 13 in ref. [50].Note that in ref. [50] 13 in ref. [50].
The data agrees over the whole coupling region, showing no evidence of a discontinuous phase transition.
convention for the electric coupling is used differing by a factor 4. Thus, the upper end of the coupling ranges is the same while our lower end goes further down to g 2 = 0. Since we also do not observe finite-size effects it supports the claim in ref. [50] that there is no discontinuous phase transition taking place.
In the second part we study fermionic observables, related to the fermionic correlations of the ground state.These are used in ref. [50] to probe a phase transition between a deconfined U(1) spin-liquid and a confined phase exhibiting antiferromagnetic order (AFM).The observable that is computed is the spin structure factor χ S (k): x,y α,β=1,2 S α β (x)S β α (y) e ik(x−y) (37) with S α β (x) = ψ † x,α ψ x,β − 1/2δ αβ γ ψ † x,γ ψ x,γ .From the spin structure factor one can compute the AFM correlation ratio defined as which quantifies the strength of AFM order (δk = (2π/L, 0) denotes the smallest momentum vector).The question addressed in ref. [50] is whether in the thermodynamic limit AFM order persists down to g 2 = 0, in other words whether the π-flux state is stable against gauge-field fluctuations.The AFM correlation ratio is computed up to lattice sizes of 16 × 16 and the crossing points between neighboring lattice sizes are extracted.
The crossing points are extrapolated to the thermodynamic limit, resulting in g 2 c,∞ = 0.15 (2).The procedure is shown in Fig. 7 which is to be compared with the Euclidean Monte Carlo study in ref. [50] where the extrapoled value is g 2 c,∞,EMC = 0.40 (5).We thus obtain  7. Benchmark for compact QED 3 coupled to N f = 2 species of dynamical fermions at half-filling: we compute the AFM correlation ratio rAFM in the variational ground state for lattice size up to 16 × 16 (left).The AFM correlation ratio is computed from the spin correlations and quantifies the strength of antiferromagnetic order.The crossing points are extracted and extrapolated to the thermodynamic limit, resulting in g 2 c,∞ = 0.15(2) (right).This is to be compared with the Euclidean Monte Carlo study in ref. [50] where also a non-zero coupling was extrapolated but at a higher value of g 2 c,∞,EMC = 0.40 (5).
qualitatively similar results in the sense that both extrapolated values are significantly larger than zero and indicate a possible phase transition but the value in our method is lower compared to ref. [50].
Another interesting quantity are the spin-spin correlations as defined in eq.(37).We compute the decay of spin correlations on a 16 × 16 lattice both in the weakcoupling region (g 2 = 0.1) and in a more strongly-coupled region (g 2 = 0.85).The result for both the full correlation function and only the connected part is shown in Fig. 8.At stronger coupling the correlation function decays to a constant value which is lower than predicted by the Heisenberg model (as to be expected since g 2 = 0.85 is still too small for a Heisenberg description).The connected correlation function decays exponentially as expected.At weak coupling the connected correlation function rather decays algebraically, as expected for a gapless spin liquid.The form of the decay is very similar to one in the Euclidean Monte Carlo study (see Fig. 4 in ref. [50]).
We can thus, at least support the claim in ref. [50] that there is indeed a deconfined phase which, however, only persists up to a smaller coupling of g 2 c,∞ = 0.15(2) in our case.One should note though that for the extrapolation of the AFM correlation ratio and also the computation of the spin structure factor is very sensitive to errors (as also mentioned in ref. [50]) so that a quantitative difference can be expected.

V. SIGN-PROBLEM AFFECTED REGIMES
In this section, we access regimes where the signproblem is present in order to demonstrate that our FIG.8. Benchmark for compact QED 3 coupled to N f = 2 species of dynamical fermions at half-filling: we compute the decay of spin correlations (both for the full correlation function and the connected correlation function) in the variatonal ground state for a 16 × 16 lattice at weak coupling (g 2 = 0.1) and at stronger coupling (g 2 = 0.85).Note that we only use odd distances in r to avoid oscillations.At strong coupling (where one expects behaviour similar to the Heisenberg model) the full correlations decay to a constant while the connected part decays exponentially.At weak coupling the decay is rather algebraically, similar to the decay shown in the Euclidean Monte Carlo study in Fig. 4 in ref. [50].
method does not suffer from the sign-poblem.Having benchmarked our ansatz for the scenario of two flavors of fermions at half-filling, i.e. zero chemical potential, it is natural to study this configuration at finite chemical potential.
We specifically want to look at a scenario that has FIG. 9. Finite potential study for compact QED 3 coupled to N f = 2 species of dynamical fermions: we compute the variational ground state energies (corrected by an overall constant µN/2) on an 8 × 8 lattice for different isospin numbers ∆N depending on the chemical potential difference µ between the two species.We study both the case of massless and massive fermions (see top row, (a,c)).By computing the crossing points between the ground state energies we can extract the phase transitions between neighboring ∆N phases (see bottom row, (b,d)).
been used in one dimension with tensor networks [33] to demonstrate overcoming the sign-problem and extend it to two dimensions.In the referenced work the authors study density-induced phase transitions due to varying flavor-dependent chemical potentials.Analogously to ref. [33], we look at the case of massless and massive fermions.
We fix the parameters in the Hamiltonian given in eq. ( 1) to the values t = 1, g mag = −1 and g 2 = 0.2, similar to the benchmarked case in the previous section.Only the staggered mass m and the chemical potentials µ 1 and µ 2 will be changed.To make this explicit we rewrite the Hamiltonian as with the conserved quantities x,2 ψ x,2 .Alternatively, one can also use the total number of fermions N = N 1 + N 2 and their imbalance (sometimes called isopsin number) ∆N = N 1 − N 2 as conserved quantities.Respectively, one defines the chemical potentials µ + = (µ 1 + µ 2 ) and µ − = (µ 1 − µ 2 ).The rest of the Hamiltonian is contained in H 0 (m) which only depends on m.
The Hamiltonian is block-diagonal and different sectors are labelled with N and ∆N .In analogy to ref. [33], we fix the total number of fermions N to the number of lattice sites and study the nature of the ground state (characterized by ∆N ) depending on µ − , the isospin chemical potential.Since the energy (up to a constant) only depends on µ − , we set µ 2 = 0 so that µ + = µ − = µ 1 ≡ µ.The ground state energy for each sector can then be written as where E 0,∆N,N (m) is the ground state w.r.t.H 0 (m) for fixed N and ∆N .We ran our simulations on an 8 × 8 lattice where we saw that finite-size effects were negligible for our purposes.To detect the phase transitions between different ∆N phases we compute the variational ground state energy for a given ∆N to determine E 0,∆N,N (m).We plot the ground state energy of every ∆N -sector subtracted by the constant given by the total number of fermions, i.e.E ∆N,N (µ, m) − N µ 2 .The crossing points between different ∆N energies give us the location of the phase transitions.
The result of that procedure for the massless case is shown in Fig. 9(a) which then allows us to plot the ∆N phase transitions, illustrated in Fig. 9(b).For the mas-sive case we choose a staggered mass of m = 1.0 and repeat the same procedure as in the massless case.The result for the ground state energy is shown in Fig. 9(c) whereas the phase transitions are shown in Fig. 9(d).
When comparing the massless and massive case it becomes clear that the phase transitions are all shifted to higher values of µ.However, the extent of this shift depends strongly on the isospin number.While phase transitions between small ∆N are severely affected by the staggered mass term (in particular the transition between ∆N = 0 and ∆N = 2 which shifted from µ = 0.5 to µ = 1.5), phase transitions for larger ∆N are relatively unaffected (e.g. the phase transition between ∆N = 6 and ∆N = 8 shifts only slightly from µ = 2.2 to µ = 2.3).This is in agreement with the results of the tensor-network study in one dimension [33].
The reason for this behaviour lies in the different changes in ground state energy E 0,∆N,N (m) for different ∆N if we go from H 0 (m = 0) to H 0 (m = 1).Qualitatively, this can be explained with the fact that for larger isospin numbers ∆N the imbalance in occupation between even and odd sites (the lattice analogue of the chiral condensate) becomes smaller and thus gets more penalized by a staggered mass term.Hence, the phase transitions shift to higher values in chemical potential.Since this effect is stronger for smaller isospin number, it mostly affects transitions between such phases.

VI. CONCLUSION
In summary, we have presented a variational, signproblem-free Monte Carlo method to study higherdimensional lattice gauge theories with dynamical fermions without truncating the gauge field Hilbert space and applied it to (2+1)-dimensional compact QED with dynamical fermions.
We benchmarked the ansatz against limiting cases of the model, against other variational methods [45] and against a Euclidean Monte Carlo study [50].To access sign-problem affected regimes we study the model at finite chemical potential, namely (in analogy to a tensornetwork study in one dimension [33]) we detect densityinduced phase transition for the case of massless and massive staggered fermions.
The variational ansatz is formulated in the gauge field basis of U (1), consisting of two parts: first, a Jastrow-type ansatz state is constructed out of the gauge field plaquette variables, thus readily gauge-invariant.It can describe the ground state of pure gauge compact QED.Secondly, a gauge-fermion part is introduced that is an infinite superposition of gauge-field dependent fermionic Gaussian states.For every gauge field configuration a fermionic Gaussian state is defined in such a way that the integral over all gauge field configurations is gaugeinvariant and still efficiently tractable, i.e. the number of variational parameters scales polynomially in system size and not exponentially.Such a construction can be achieved by making the variational parameters of a fermionic Gaussian state gauge-field dependent and use a parametrization based on the eigendecomposition of the gauge-matter Hamiltonian.This requires exact diagonalization at every measurement step in the sampling algorithm but since its scaling is O(N 3 ) in system size, the method is efficient and we can reach large lattice sizes.
In the future it would be interesting to also study other higher-dimensional lattice gauge theories such as three-dimensional or non-Abelian lattice gauge theories.The former would not require any change in the ansatz while the non-Abelian nature of the gauge group only requires changes in the pure gauge part of the ansatz.The fermionic part could stay the same since the eigendecomposition of the gauge-matter Hamiltonian can still be carried out efficiently.The pure gauge part would need to be changed since the plaquette operators on which the Jastrow wavefunction is based would not be gauge-invariant anymore but this could be remedied by using traces of plaquette operators and potentially other closed loops.
Also, in light of more and more powerful quantum devices which require for the simulation of lattice gauge theories some kind of truncation in the gauge field Hilbert space, the presented method could help in studying and thus controlling the errors caused by such truncations.
For the case of N f = 2 fermion flavors (and potentially even more flavors) one can choose a similar block-diagonal structure where one now blocks the single-particle eigenstates of both flavors together, i.e. ψ † 1,i+ (θ) |0 , ψ † 1,i− (θ) |0 , ψ † 2,i+ (θ) |0 , ψ † 2,i− (θ) |0 .The individual blocks ξ| i are then 4 × 4-matrices (or more generally 2N f × 2N f -matrices).The variational parametrization of these blocks is kept general as this allows to control certain properties between the two species, e.g. the imbalance ∆N between the two species as was used for the study of sign-problem affected regimes in section V.The block-diagonal structure allows even for multiple fermion flavors to compute the local energy H loc (θ), in particular H E,loc (θ), (see section III B 1) with a computational cost of only O(N 3 ).

2 FIG. 4 .
FIG. 4. The flux θp per plaquette is shown for the variational ground state at g 2 = 0.0, t = 1 and gmag = −1 for a 12 × 12 lattice.The average deviation of θp from π is on the order of 10 −4 .The global loops θ1 and θ2 (winding around the axes of the lattice) also acquire a π-flux with a similar deviation as the plaquette fluxes.
FIG.5.Benchmark for pure gauge compact QED 3 of our ansatz (denoted by FG) against the variational method in ref.[45] based on complex periodic Gaussian states (denoted by CPG): we compare the ground state energy E0 for an 8 × 8 lattice over the whole coupling region and compute the relative error (see inset).
FIG.7.Benchmark for compact QED 3 coupled to N f = 2 species of dynamical fermions at half-filling: we compute the AFM correlation ratio rAFM in the variational ground state for lattice size up to 16 × 16 (left).The AFM correlation ratio is computed from the spin correlations and quantifies the strength of antiferromagnetic order.The crossing points are extracted and extrapolated to the thermodynamic limit, resulting in g 2 c,∞ = 0.15(2) (right).This is to be compared with the Euclidean Monte Carlo study in ref.[50] where also a non-zero coupling was extrapolated but at a higher value of g 2 c,∞,EMC = 0.40(5).

1 S
(r)S(0) conn.corr.(g 2 = 0.10) full corr.(g 2 = 0.10) conn.corr.(g 2 = 0.85) full corr.(g 2 = 0.85) Benchmark for compact QED 3 to N f = 2 species of dynamical fermions at half-filling: the shown averaged flux energy per plaquette cos(θp) in the variational ground state is to be compared with results obtained in an Euclidean Monte Carlo study shown in Fig.