Resource-Efficient Quantum Simulation of Lattice Gauge Theories in Arbitrary Dimensions: Solving for Gauss' Law and Fermion Elimination

Quantum simulation of Lattice Gauge Theories has been proposed and used as a method to overcome theoretical difficulties in dealing with the non-perturbative nature of such models. In this work we focus on two important bottlenecks that make developing such simulators hard: one is the difficulty of simulating fermionic degrees of freedom, and the other is the redundancy of the Hilbert space, which leads to a waste of experimental resources and the need to impose and monitor the local symmetry constraints of gauge theories. This has previously been tackled in one dimensional settings, using non-local methods. Here we show an alternative procedure for dealing with these problems, which removes the matter and the Hilbert space redundancy, and is valid for higher space dimensions. We demonstrate it for a $\mathbb{Z}_2$ lattice gauge theory and implement it experimentally via the IBMQ cloud quantum computing platform.


I. INTRODUCTION
Gauge theories, describing the fundamental interactions among the constituents of matter, pose a serious challenge. Many are non-perturbative, at least for some energy scales; e.g., Quantum Chromodynamics (QCD), the theory of the strong nuclear force, is asymptotically free in high energies [1, 2], but non-perturbative at low energies. Perturbative techniques fail to describe this regime which exhibits very important physical phenomena, such as quark confinement [3] which is responsible for the hadronic structure. This strongly coupled physics has been successfully addressed (see e.g. [4]) by applying Monte-Carlo methods to lattice gauge theories (LGTs) [3,5,6] -lattice formulations of gauge theories. However, these methods cannot directly describe real-time dynamics (being based on Euclidean time) or the physics of fermions with finite chemical potentials, due to the sign problem [7]. Thus, quantum simulation [8], where a hard-to-solve quantum problem is mapped to a highly controllable quantum device which can be studied experimentally, would be useful in this case.
First, the matter is usually fermionic, and the gauge field is not. This entails combining fermionic and non-fermionic ingredients in the simulator, which often leads to choosing ultracold atomic systems for simulations [10]. In a single space dimension (d = 1) this can be overcome using the Jordan-Wigner map [24] which replaces fermions by spins -but introduces non-locality. Second, LGTs are highly constrained: the local symmetry introduces conservation laws (Gauss' law) on every site, giving rise to a redundancy in the Hilbert space, which requires the simulation of unnecessary degrees of freedom, wasting costly resources. Finally, in d > 1, the LGT Hamiltonian introduces the nontrivial four-body plaquette interaction [5], which is not possessed naturally by the common quantum devices.
The latter issue may be dealt with in several ways, and in particular using a Trotterized approach [25,26]: instead of mapping the Hamiltonian of the simulated model to that of the simulator (which is known as analogue quantum simulation), one approximates the time evolution operator exp (−iHt) = exp −i i H i t by a sequence of short time unitaries: exp (−iϵH i ), for N = t/ϵ large enough such that (1) By implementing each H i individually, complicated interactions may be composed out of two-body unitaries, possibly using auxiliary ingredients. This is useful, in particular (but not only), for the four-body plaquette interactions included in LGTs [27][28][29][30][31][32].
Here we address the first two issues by using a reformulation of lattice gauge theories which uses the local constraints to eliminate the fermionic matter [33,34]. Here we address the first two issues by using a reformulation of lattice gauge theories which uses the local constraints to eliminate the fermionic matter [33,34]. We construct a quantum simulation algorithm based on it, show specifically how to apply it to Z2 LGT to build a working quantum simulation, and demonstrate an experimental implementation of it. Our protocol yields a local model that does not include matter, but is still equivalent to the original LGT (with fermions); as a result, not only do we not need to simulate any fermions directly, but there are also no local constraints to impose and maintain and no redundancy in the Hilbert space. We thus obtain a much simpler simulation scheme that is valid also for d > 1.
For the sake of completeness, we would like to mention that other methods, using different types of tools, are also used in order to address the three issues mentioned above. These include, for example, the loop-string-hadron formalism, which formulates lattice gauge theories in terms of explicitly gauge arXiv:2206.00685v3 [quant-ph] 8 Aug 2023 invariant degrees of freedom, allowing one to remove the constraints [35,36] (see [37] for a comparative study of this and other methods for the quantum simulation of S U(2) models in a single space dimension). Another approach is the use of dual formulations, which allow, at least in the Abelian case, to switch to magnetic degrees of freedom, which are free of the Gauss law constraints and have no plaquette interactions, but do not directly address the issue of fermionic matter [38][39][40][41][42][43].
We would also like to mention, that after the completion of the first version of this article, we became aware of a parallel work on simulating Z 2 lattice gauge theories using the same methods [44]. The analysis performed in the two works may be seen as complementary.
The article is organized as follows. We begin by reviewing the basics of Hamiltonian LGTs (section II A), focusing on Z 2 as the simplest case, which already shows the relevant features (redundancy of the Hilbert space and fermionic matter). In section II B, we review the conventional non-local, 1 + 1d techniques to deal with these issues [17,[45][46][47][48] and apply them to Z 2 . The main step involves a unitary transformation that we introduce and denote as U (0) . Section III introduces our local procedure and applies it for Z 2 in d = 1. The procedure involves two unitary steps that we introduce and denote as U (1) and U (2) . We proceed to presenting a few experimental demonstrations (section IV) of the method implemented on IBMQ devices. Next, we generalize our procedure to d = 2 (section V), and present an experimental implementation of a quasi-two-dimensional system (section VI) which is the simplest system for which the standard method of U (0) fails. Finally, in section VII we use numerical simulations to estimate the near-term experimental feasibility of using our method for more advanced applications than those presented in section IV.

A. Hamiltonian lattice gauge theories
Hamiltonian LGTs [5] are defined on d dimensional spatial lattices Z d (square/cubic by default, though formulations in other geometries exist).
LGTs include two types of fields: the matter, mostly (but not necessarily) fermionic, associated with the lattice sites and described by the fermionic Fock space H m , and the gauge field, associated with the lattice links and described by the Hilbert space H g (see Fig. 1). The gauge group G is a compact Lie or a finite group that generates gauge transformations: local unitaries Θ g (x) under which the physically relevant states and operators are invariant. These are parametrized by group elements g ∈ G, and are associated with the sites x ∈ Z d . Each Θ g (x) acts locally on x and the links ℓ ∋ x around it (starting or ending at x, see Fig. 1), transforming only those degrees of freedom in a way that is parametrized by g. A gauge invariant operator O satisfies  (7)).
non-Abelian case, gauge transformations can mix the elements of state multiplets [49]). Let us focus on the case G = Z 2 . Each site can host a single fermion, annihilated by ψ (x). Each link hosts a two dimensional Hilbert space. Z 2 contains a single nontrivial group element; thus the possible gauge transformations are where acting on the links ℓ that are connected to the site x, and N (x) = ψ † (x) ψ (x) is the number operator at x. The gauge invariant operators are Z operators, products of X operators along closed loops, and functions thereof; but also (functions of) the so-called mesonic strings, which are operators of the form where C is any path connecting the sites x,y. Note that the mesonic strings include trivially the number operator N (x). A conventional Hamiltonian choice [5] takes the form where H E , the electric energy, is a sum of local gauge (electric) field terms on all the links ℓ, the magnetic energy H B , is a four-body interaction of the links around each plaquette (unit-square), and H GM is the interaction with the matter that involves hopping of fermions to neighbouring sites, while changing the state of the field on the intermediate link. H m is a the mass term, which we choose to be staggered [50], with generalizations to HEP-like LGTs in mind. While the Hamiltonian was originally formulated by Kogut and Susskind for continuous groups [5], following Wilson's Lagrangian formalism [3], it is possible to extend it to finite groups, and in particular to Z N . Here, we follow the formulation of [51], where the pure-gauge parts of the Hamiltonian are constructed such that they will give rise to the conventional U(1) formulation in the large N limit. For the Z 2 case the Hamiltonian terms can be written as: where the indices 1-4 label the four different links that form a given plaquette p (see Fig. 1), where e i is a lattice vector in direction i, and X (x, i) acts on the link emanating from x in direction i; and where the alternating sign is a result of staggering, such that for the odd sites the existence of a fermion (N (x) = 1) can be interpreted as the vacuum (Dirac-sea) state, and the absence of a fermion (N (x) = 0) can be interpreted as an anti-particle. The even sites follow the opposite and more intuitive convention where N (x) = 0 represents the empty state and N (x) = 1 represents a particle [50]. Gauge invariant states |ψ⟩ satisfy the local Gauss' law constraints, that for Z 2 can be written as: Since H is gauge invariant, the eigenvalues q (x) = 0, 1 are constants of motion, splitting the Hilbert space H into dynamically disconnected superselection sectors, H ({q (x)}). The physical Hilbert space satisfies In a model with staggered mass as in Eq. (9), the sector defined by e iπq(x) = (−1) x 1 +...+x d is often considered as the simplest sector in the sense that it includes the "Dirac-sea" state in which only odd sites are populated (no particles and no antiparticles).
B. Removing the redundancy in the standard approach: Eliminating the gauge fields Removing the redundancy means solving Gauss' laws (10): for every lattice site x. In the standard approach, one solves it for the gauge fields: given N (x), we have to solve for Z on each link. However, there are several Z configurations satisfying Eq. (12), unless d = 1 (Fig. 2(a)). In this case, we can label both the sites and the links by n, and the constraints simplify to Z n Z n−1 |ψ⟩ = e iπ(N n +q n ) |ψ⟩ (13) which, for open boundary conditions (0 ≤ n ≤ L − 1, for an even number of sites L) is easily solved by the non-local expression: This motivates the unitary degauging transformation: We denote transformed states and operators as This transformation enforces the solution (14) in a given sector. The transformed parts of the Hamiltonian are the interaction: and the electric term, which becomes non-local: In d = 1 there is no H B (no plaquettes), and the mass term is unaffected: Note that H (0) , Z n = 0 and Z n ψ (0) = ψ (0) , ∀n. Thus the redundancy is completely removed, the gauge fields are in a product state with the matter, and only the latter has to be treated in any quantum simulation scheme. The explicit presence of fermions is still a potential problem, restricting the choice of the simulating platform, but in d = 1 we can use the Jordan-Wigner transform [24] to represent the fermions with qubits: which results in an L qubit system, residing on the sites ( Fig.  2(b)). The Hilbert space dimension has decreased exponentially from 2 2L−1 to 2 L , with no local constraints left. This reduction is demonstrated in Figs. 3(a,b), by comparing the spectra of H and H (0) . As H (0) is in a specific sector it contains less levels, but it describes the same physics as H in the chosen sector. For simplicity, we will focus from now on the sector where e iπq n = (−1) n , in which the Hamiltonian is (up to a constant): The dynamics of H (0) m can be simulated using local qubit rotations, and that of H (0) GM by simple two-qubit gates on neighbouring sites. H (0) E , however, involves highly non-local manybody interactions, and thus it is more complicated and requires Trotterization, either to a strictly digital simulation or an analogue-digital one. In the latter, the evolution with respect to H (0) m and H (0) GM can be implemented using analogue techniques (which is possible on some platforms). In any case, the local parts can be simulated with an O(1) run-time.
To implement e −iϵH (0) E , we use two types of unitaries: (i) single qubit rotations, V n = exp −iϵh (−1) n(n+3)/2 σ z n ; (ii) CNOT gates between neighbouring qubits, U n = |↑⟩ ⟨↑| n + |↓⟩ ⟨↓| n ⊗ σ x n+1 . Since U n σ z n+1 U n = σ z n σ z n+1 , we get that Due to the non-locality, the length of this operation scales linearly in the system size L, and we conclude that quantum simulation of the dynamics using this method has an O(L) runtime per Trotter step.
While we demonstrated it for Z 2 , this way of integrating the gauge field out is valid for arbitrary gauge groups in d = 1 [34,[45][46][47], and can be used for quantum simulation [17,48]. This method's drawbacks are that it is restricted to d = 1, introduces non-locality and its Trotter step run-time depends on L. non-locality can arise in different ways, for example for Lie groups -see, e.g. [17], where the transformed electric Hamiltonian includes only two-body terms, but arbitrarily far; in that case, a successful experimental realization was possible using the long-range interactions of the simulating platform used (trapped ions).

III. ELIMINATING THE MATTER
In our approach we solve the constraints as equations for the matter. This significantly simplifies the problem, since in this view these equations are explicitly solved: In the Z 2 case, knowing the Z configuration gives rise immediately to a unique and local solution for N (x), and similar results are valid for other groups as well. Such a solution is quite straightforward for bosonic, Higgs-like matter (unitary gauge fixing) [61]. When we deal with fermions, things have to be done rather more carefully, but it is possible nevertheless. There are two steps to our procedure: first (section III A) we transform the fermions into hard-core bosons, and then (section III B) we solve Gauss' law for the matter and remove the redundancy, eliminating altogether the need to simulate the matter. In section III C we provide details on how to implement the quantum simulation (transformed) Hamiltonian using fully digital or hybrid techniques, as well as how to measure the gauge invariant observables.

A. From fermions to hard-core bosons
Following the procedure of Ref. [33] we can replace the fermionic matter of any LGT whose gauge group contains Z 2 as a normal subgroup by hard-core bosonic matter. After applying a unitary procedure U (1) which preserves the physics of the original states |ψ⟩, as given in [33], one ends up with equivalent states, of a model in which each fermionic mode ψ is replaced by a spin, or a hard-core boson; thanks to the local constraints (Gauss' law), the gauge field absorbs the statistics, leaving us only with non-fermionic matter fields and modifying slightly the way that the gauge fields appear in the Hamiltonian, to account for the statistics. Since this is done via local unitaries which exploit the local constraints, we are left with a local theory nevertheless [33]. While, as shown in [33], the method is valid for any space dimension we will focus here on d = 1 for the clarity of the presentation, and for comparison with the previous method (we treat the d = 2 case in section V). On the other hand, we will switch from now on to periodic boundary conditions, which are simpler to deal with and, unlike when the gauge field is removed, are valid here. For Z 2 , when performing the procedure of [33] and replacing the fermions by hard-core bosons, we get where σ ± n are spin raising and lowering operators for the hardcore bosonic mode at the site n.
Transforming the original Gauss' law (given by Eq. (13)) with the substitution N n → σ + n σ − n , we obtain its hard-core bosonic form, where S n = Z n−1 Z n , and ε n ≡ e iπq n = ±1 such that each choice of signs {ε n } L−1 n=0 defines a superselection sector. The Hamiltonian H (1) is free of fermions, but is subject to local constraints. The redundancy problem is unsolved, but one can still formulate a Trotterized quantum simulation of H (1) , using 2L qubits (2L − 1 for open boundaries). The usual set of gates and tools allows us to formulate it quite easily, but we still need to use a redundant Hilbert space and make sure that the constraints are satisfied, either by monitoring them directly [57,62,63] or making sure that each Trotter step is gauge invariant [27,31]. In the recent work [57], the d = 1 Z 2 was simulated without integrating out any degree of freedom (the fermions were taken care of by a Jordan-Wigner transform instead, and thus extending it to higher dimensions might be very challenging and non-local. Such ideas have been studied recently, e.g. in [64] and references therein).

B. Solving Gauss' law for the matter
Here, we shall proceed to a complete elimination of the matter, in a procedure similar to that of Ref. [34], using the fact that Gauss' law provides us with a one-to-one map between the values of S n and those of σ z n . An alternative approach for eliminating the matter, valid for Z 2 only but giving rise to similar results, was studied in [65,66]; similar methods for eliminating matter in d = 1 abelian systems were discussed in [67][68][69]. Our procedure is valid for other gauge groups (including non-Abelian ones) and higher dimensions as well. Originally, it was given for U(N) groups, and we shall now adapt it to the case of Z 2 , which was not explicitly included in [34].
First, we define projectors onto S n eigenstates where ε n = ±1, depending on the static charge sector. and rewrite Gauss' law as valid for any choice of signs {ε n }. Define, on each site n, a controlled local unitary which decouples the matter: Suppose we want all the matter spins to point down. If S n = ε n , we do nothing, and if S n = −ε n we invert it (refer to Eq. (25)). The controlled unitary that performs this operation is: and since [U n , U m ] = 0, we can safely define This is the second unitary step in our procedure, and we denote transformed states and operators as and Note that the operators U (2) n depend on the projection operators P ± n which depend on the static charges. Hence, our transformation is valid for a given sector on the Hilbert space, or in other words, constructed to fit a the sector of interest. Thanks to the superselection of static charges, there is no point in discussing more than a single sector, and this is the point where we make an explicit choice of the sector, discarding all other sectors henceforth.
By construction, the matter qubits are completely decoupled in the transformed state, as the transformed Gauss' law (apply U (2) to Eq. (27)) is: -transformed physical states are ones in which all matter qubits are in the σ z n = −1 state. In other words, we started with a state with gauge fields and matter, satisfying the local Gauss' law constraints, and ended up with a state where the gauge fields and matter are decoupled, and the local constraints are satisfied by the matter degrees of freedom alone. Originally, the Hilbert space was divided into dynamically disconnected sectors given by Gauss' law, and now the sectors are of the decoupled matter alone ( H (2) , σ z n = 0 ∀n). For this reason, in the beginning, while being constrained, we could not simply discard the matter degrees of freedom, now it possible to do so thanks to the decoupling.
We can thus restrict ourselves to the sector where all the matter spins point down. They are not affected by the dynamics, and hence they do not have to be simulated. Formally, if we define by |out⟩ ∈ H (2) m the matter-state for which σ z n = −1 ∀n, our relevant quantum simulation Hamiltonian will bẽ H (2) acts only on field (link) qubits states and describes the same physics as the original H in a specific chosen charge-sector defined by the choice of {ε n }. We thus arrive at a theory in a much smaller Hilbert space, but with no local constraints, containing only the relevant part of the spectrum (as can be seen in Fig. 3). This is the result of combining a unitary transformation (preserving the spectrum) and a projection (which keeps only the relevant part of it). By simply plugging different static charges to the definitions of the projectors and the transformation, one can obtain a similar result for any other sector. Importantly, the choice of matter sector (on our case -all matter spins pointing down) is not completely orthogonal to the choice of charge-sector {ε n }, and one has to check for consistency with the global charge symmetry, where q = n q n . Since the sign of the right hand side is determined by the static charge sector and the left hand side by the fermionic parity sector, the two choices have to be made such that Eq. (34) is fulfilled. In practice this means that for charge sectors with an odd q, one would have to use a slightly different matter sector instead of the one we use here (for example -the first matter spin in the chain points up, and all the others point down), and the decoupling operation U n would have to be changed accordingly. Applying this procedure to H (1) (Eq. (24)), one finds that the electric term remains unchanged, the mass term becomes the local two-body interaction: and the interaction term takes the form This procedure (applying U (2) to H (1) and projecting on |out⟩) is described in more detail in Appendix A. At this point we focus an the specific charge-sector with ε n = (−1) n (chosen to include the "Dirac-sea" state). This is consistent with our matter-sector choice only when L is an integer multiple of 4, so from here on we restrict ourselves to this case. Note that the procedure can be easily altered to fit the other even L case instead (by choosing a different mattersector, and changing U n accordingly as explained above). As a result, the mass term simplifies, butH (2) GM still has an alternating sign (−ε n ) = (−1) n+1 . This is not a problem, but for the sake of elegance we make one extra step, using to finally obtain: and we have re-defined the interaction and electric parts such that the former includes only three-qubit interactions and the latter has all the single qubit terms (including those that came from the original interaction part).
To change to open boundary conditions one can use almost the same expressions, but remember to sum over the sites 0 ≤ n ≤ L − 1 forĤ m , and over the links 0 ≤ n ≤ L − 2 for H E andĤ GM . Then one has to make the substitution Z −1 = Z L−1 = 1 which can be thought of as placing two additional links at the boundaries, with fixed field values. The result is the addition of boundary terms that are simpler than the bulk terms (single-qubit instead of two-qubit terms, and two-qubit instead of three-qubit terms).
In both cases, we now have an L−1 link-qubits Hamiltonian (Fig. 2(c)), acting on states ψ = V |ψ (2) ⟩ without constraints, global or local: again, an exponential reduction of the Hilbert space (see Fig. 3(a,c) for a comparison of the spectra of H andĤ). However, in contrast to the standard method of H (0) (section II B), we now have a local Hamiltonian, and the procedure is generalizable to higher dimensions (see section V).

C. From Hamiltonian to quantum simulation
Time evolution with respect toĤ is readily implemented using Trotterization. Ω E (ϵ) ≡ e −iϵĤ E (where ϵ is the length of a Trotter step) can be implemented in an analogue fashion (that is, as a whole) during the Trotter step; on the other hand, in a more digital approach, it can be decomposed into a product of local, commuting single qubit rotations, where r = h 2 + J 2 /4, and cos θ = −h/r and sin θ = −J/2r define the axis of rotation in the ZY plane. It is very likely that any simulating platform will be able to run all these gates in parallel, and even if not, it should be possible to do it in a finite number of steps where several qubits are rotated in parallel. Thus, for all practical purposes one can assume that Ω E is implemented in an analogue way. A similar argument holds for Ω m (ϵ) ≡ e −iϵĤ m . Here, instead of local terms we have two-body ZZ interactions of nearest neighbours (and local rotations for the ends of the open system) which mutually commute, and may be implemented either together (analogically -note thatĤ m is nothing but a simple Ising Hamiltonian) or sequentially with a small number of steps (since some of the constituent operations can be run in parallel, depending on the simulating platform).
Finally, Ω GM (ϵ) ≡ e −iϵĤ GM would be more challenging for most simulation platforms, since it involves three-body interactions which are not natural for them. To implement it, we use the conventional controlled-Z gate: which obeys Defining U CZ = n U CZ n , it follows from Eq. (44) that where U Y (ϵ) = exp iϵ J n Y n /2 ≡ exp (−iϵH Y ) which, again, can be run either in parallel or sequentially, depending on technological constraints of the simulating platform. In either case, it can be done with a finite number of steps, independent of L, and we conclude that the entire algorithm runs in O(1) time.
The only remaining task is to choose the order of the three unitaries out of which a Trotter step is built. A Trotter error analysis, which is given in Appendix B, shows that the optimal ordering is which can be applied as a recipe for a fully digital quantum simulation of the model. The exponential form of the CZ operation can be used to construct a hybrid analogue-digital simulation: First, use Eq.
Since H Z andĤ m not only commute, but also have a very similar functional form, we can definê which is a simple Ising Hamiltonian with a longitudinal field. Then we can obtain our single Trotter step using a sequence in which we switch on and off four analogue Hamiltonians: After simulating time evolution (either in the hybrid or in the fully digital way), we have to be able to measure observables from the original model. The relevant local observables are the electric field on the links, and the number operator N n = ψ † n ψ n on the sites. To measure these, we first have to check how they transform under our procedure: first with U (1) , then with U (2) , and finally projecting the matter state onto |out⟩ and rotating with V (though in these cases V has no effect). It is easily verified that under this procedure E n and N n transform to: The field E n is unchanged and can therefore be obtained trivially from measuring the qubits in the computational basis, while for N n we have to measure the product S n = Z n−1 Z n , which is the parity of neighbouring qubits. The non-local gauge invariant observables are the mesonic strings, defined in Eq. (4). Measuring those is also possible within this scheme, but it is somewhat more involved and we show how to do it in Appendix C. This concludes our construction for the one-dimensional case. Such a simulator would be useful for a broad range of tasks, e.g. adiabatic ground state preparation, or studying quenches, some of which are exemplified in the following section.

IV. EXPERIMENTAL IMPLEMENTATION
We implemented a proof-of-concept version of this quantum simulation proposal via the IBMQ platform. For that we focus on the 1 + 1d case with m = 0 and open boundary conditions, in the ε n = (−1) n sector. This means that we have to implement the L − 1 qubits Hamiltonian:  where the last two terms are boundary terms. The hybrid analogue-digital approach (Eq. (49)) might possibly be implemented on those IBMQ devices that allow for direct pulse control, but this is beyond the scope of this work. Instead we follow the Trotterization procedure for a fully digital simulation, which can be summarized by Eq. (42) and (45) The operation U CZ (controlled-Z on all pairs of neighbouring qubits in the chain) that appears twice in Ω GM (ϵ) has to be implemented in two steps (one for the even pairs and another for the odd pairs). This means that each Trotter step can be implemented with 4 two-qubit gate steps, and 2 singlequbit rotation steps. We emphasize again that these numbers do not depend on L. Converting from CZ gates and general single-qubit rotations to the native gates of the IBMQ devices (CNOT, X, √ X and virtual Z gates) costs in additional 4 single-qubit steps.
Typical IBMQ qubits have coherence times on the order of 100 microseconds, and native single-qubit gates can be implemented within 35ns. Two-qubit (CNOT) gates however, are implemented with via the cross-resonance approach [70,71] and typically take between 300-500ns each. Assuming we want the computation to complete within ∼ 10% of the coherence time, this restricts us to about 5 Trotter steps in total (about 20 native two-qubit steps and 30 native single qubit steps where each native step acts on the entire chain). This poses a limitation on the possible computations. For example: when implementing adiabatic ground-state preparation, the adiabaticity condition cannot be fulfilled for some regions in parameter space, resulting in poor fidelities. Nevertheless, it is important to remember that faster or more coherent hardware does exist, and state-of-the-art technology already allows for an order of magnitude improvement in the coherent Trotter depth. The rather strict requirement of completing the experiment within 10% of the coherence time is an empirically (and numerically) verified heuristic that seemed to optimize the Trotter error against decoherence errors in most of our experiments. However, it has been shown that error mitigation techniques like ZNE (which we did not implement here) allow for meaningful evaluation of observables even when a larger degree of decoherence noise is allowed in the experiment [72].
One of the most significant advantages of quantum simulation is the possibility to simulate time evolution. Importantly, the limitation of 5 trotter steps does not translate to a limitation on the temporal resolution, since one can directly control the size of each step (which translates to an angle of rotation in a single-qubit gate). Practically this means that we have to choose the number (between 1 and 6 in this case) and the length (between 0.4/J and 0.5/J) of the Trotter steps to fit each desired simulated evolution time. For our demonstration we initialize the L = 4 chain with an excitation in one of the qubits: this corresponds to an excitation of the field on the relevant link, as well as a change in the sites connected to it to accommodate the original gauge constraints. Then we evolve it in time and measure the local observables (E n on the links and N n on the sites) as a function of the evolution time. The measurement is averaged over 20000 to 30000 shots such that the readout error is insignificant. We observe (Fig. 4) that for small values of h/J the initial excitation diffuses to the neighboring sites and links, while for large h/J it remains confined. With only 4 sites, we cannot claim to having observed a phase transition, however this is still a non-trivial physical feature of the model that our quantum simulation captures using only 3 qubits and a few tens of noisy gates.
Quantum simulation can also be used to investigate nontrivial ground-states via adiabatic ground state preparation. For example, since the ground state of the J = 0 Hamiltonian is trivial (all qubits are at |0⟩, which corresponds to the "Dirac-sea" state of the original model), by running a time evolution experiment while increasing J from zero with each step, we can measure the ground state for a finite J.
Motivated by recent work on scaling phenomena near the h = 0 transition [73], we chose to implement the opposite (increasing h adiabatically) for the purpose of our proof-ofconcept demonstration. Initializing the h = 0 ground state is not as trivial as the J = 0 ground state, but there is a simple shallow circuit that initializes the ground state ofĤ GM (that is, only the term that is interacting for the qubits, rather then the interaction term of the original model). This circuit is straightforward to derive based on Eq. (44).
After this initialization we proceed by adiabatically increasing the non-interactive terms simultaneously, and arrive at the desired finite h ground state. This scheme was implemented for L = 4 sites and the results are summarized in Fig. 5, showing good agreement with the exact solution and with a noisy numerical simulation, implemented on Python via the Qiskit-Aer package. For transparency, we used a custom noise model that includes only energy-relaxation and dephasing channels, with T 1 , T 2 for each qubit and duration for each gate as reported by IBMQ.
Thus, we can probe the ground states of both the small h and the small J regimes. Intermediate regimes are more challenging on the IBMQ devices due to the aforementioned limitation on the total number of Trotter steps, but we show numerically (Fig. 8) that reasonable fidelities can be expected with current technology. This is discussed further in section VII.

V. GENERALIZATION TO TWO SPATIAL DIMENSIONS
As we showed in section II B, the traditional methods that treat the Hilbert space redundancy and the problem of simulating fermions cannot be extended beyond d = 1. The reason for that is that the Jordan Wigner transformation assumes an order over the sites, which in d > 1 would have to be defined in an arbitrary way that is highly non-local and does not respect the lattice geometry. This is possible to in principle but extremely impractical. Even worse -the construction of U (0) relies an the existence of a well-defined solution (Eq. (14)) of the constraints for the gauge-field, which is not available in d > 1. In contrast, our procedure is completely local, and relies on a unique solution of the constraints for the matter, which is available in any dimension. We demonstrate it here for Z 2 with d = 2, in the charge sector defined by a choice of signs First, consider the hard-core bosonic formulation of the model. Applying the procedure of Ref. [33] to the Hamiltonian (5) at d = 2, we get (since the terms get rather complicated in terms of coordinates and directions, we show it graphically): Gauss' law is very similar to that of Eq. (25): with S (x) as defined in Eq. (3), completely analogous to S n in d = 1. Here we see again that when treating Gauss' law as an equation for the matter (σ z (x)) rather then for the field, it is explicitly solved, and extending to d > 1 does not change that. Therefore we can similarly define and rewrite Gauss' law as The local controlled unitaries are defined the same way: and since they all commute we can safely define U (2) = x U (x), from which the decoupling of matter follows, in the form of the new constraints From this, we can obtain the d = 2 Hamiltonian in a similar manner; it will involve local few-body interactions (since H (1) is local, and U (2) is local) which can be implemented using the same digital or digital-analogue tools. The locality guarantees that the Trotter steps can be concluded with an O(1) run-time, as in the d = 1 case. To get an idea of the result, we again focus on the simple charge sector ε (x) = (−1) x 1 +x 2 , and restrict ourselves to J ∈ R, which allows for some simplification in the resulting expressions: Here, too, we can remove the staggering with a unitary V which acts with Z on all the links emanating from sites for which x 1 + x 2 is even, and get the simulation Hamiltonian: which has the exact same terms as in Eq. (64), but without the alternating signs in front of the J/2 terms. This is, as expected, a Hamiltonian involving local qubit interactions, which can indeed be simulated using the usual quantum simulation techniques, such as those used for d = 1 in section III C. Moreover, one can repeat the entire procedure in the same way for d > 2: H (1) will have a slightly different form, but nevertheless local, and this will be the only significant change. All the arguments and techniques from section III C remain valid, and one is able to construct Trotter steps with O(1) runtime.

VI. EXPERIMENTAL IMPLEMENTATION OF A QUASI TWO-DIMENSIONAL SYSTEM
In order to implement the 2+1d version of our proposal (Eq. (65)) the qubits have to be organized on a square lattice. As this is not the case for any IBMQ machine, we implemented a quasi two dimensional version of the model with 4 sites as depicted in Fig. 6, where the middle site is treated as an "odd" site for the purposes of staggering and choosing a superselection sector. This means that the Gauss' laws on the four sites are Z n |ψ⟩ = e iπN n |ψ⟩ , for n = 0, 1, 2, FIG. 6. The quasi two-dimensional system with four sites, and the indexing convention for the sites (in purple) and for the links (in green). The middle site (n = 3) is considered an "odd" site in our chosen sector, which implies that (a) the J = 0 ground-state is the one where N 3 = 1 and all other N n and E n equal zero (as indicated by the purple highlighting of the middle node). (b) The initial state of the time evolution experiment (Fig. 7), with excitation in qubit (link) 1, is the one where N 1 = E 1 = 1 and all other N n and E n equal zero (as indicated by the highlighting of node 1 and link 1).
which implies that at J = 0 the ground state is the one where E n = 0 and N n = 0 for n = 0, 1, 2, and N 3 = 1). This toymodel is the simplest system where the standard approaches for eliminating the fermions fail due to the dimensionality and the connectivity. Assuming m = 0 as in section IV and following the matter elimination procedure, we find that the electric term H E does not change, and the interaction term H GM becomes: Plotted is the excitation with respect to the J = 0 ground state ( Fig. 6(a) ): that is, for the links we plot the field ⟨E n ⟩ = 1 2 (1 − Z n ) , for sites n = 0, 1, 2 we plot ⟨N n ⟩ and for the middle site (n = 3) we plot ⟨1 − N n ⟩ that can be thought of as the number of anti-particles. The initial state is the one where qubit 1 is excited, which corresponds to the original-model state shown in Fig. 6(b). For large h/J the initial state is more robust to the dynamics. This dynamics can be Trotterized with 8 two-qubit (CZ) gates per Trotter step and about 10 single-qubit gates (the details are in Appendix D), which means that on IBMQ machines we can preform only 2 or 3 Trotter steps within 10% of the coherence time. Unfortunately this is not enough for adiabatic ground state preparation with acceptable fidelities, so we focus on time evolution with an initial excitation (similar to Fig. 4), that allows us to observe qualitative features of the original model. In this experiment we begin by exciting qubit 1, which corresponds (in the chosen sector) to an initial state with E 1 = N 1 = 1, and E 0 = E 2 = N 0 = N 2 = N 3 = 0 (see Fig. 6). The resulting time evolution (Fig. 7) is similar to the one-dimensional case in the sense that again we observe different behavior for different values of h/J given the same initial excitation, whose robustness to the dynamics may be qualitatively related to confinement or deconfinement.

VII. DISCUSSION AND SUMMARY
In order to assess the feasibility of our method, we discuss the hardware requirements for more advanced implementations. First, in order to simulate the truly 2 + 1d version of the model, the qubits have to be arranged on a square lattice. As such devices already exist (e.g. the famous Google Sycamore [74]), this particular technical problem can be considered solved. Next, we consider what are the requirements on the quality of the qubits for extending beyond the proofof-concept demonstrations shown in this work. For that we numerically simulate the adiabatic ground-state preparation experiment described above (section IV) -with a trivial initialization of the J = 0 ground-state and an adiabatic increase of J through the Trotter steps up to some finite value.
We do this using different values for the qubits' coherence time and the CZ gate-duration, and record the fidelity of the final state relative to the exact ground state obtained by diagonalization. Taking J = h = 1 and different values of L, we can conclude ( Fig. 8(a)) that while going beyond L = 4 is challenging for the IBMQ devices, it should be possible with slightly longer coherence times or faster gates that are achievable on other platforms [75]. While a physical interpretation as a Z 2 LGT is only valid when L is an integer multiple of 4 (see section III B), evolving the simulation Hamiltonian (53) for other values of L is still a valid approach to assessing how well the method scales with the system size. Similarly, Fig.  8(b) shows the hardware requirements for probing the L = 4 model with intermediate values of J/h (recall that it is possible to start the adiabatic ground-state preparation procedure from either the h = 0 or J = 0 direction). The jumps in the plots are an artifact of the choice of a different number of Trotter steps for each data point in an attempt to optimize Trotter errors against decoherence errors. In principle, the total simulated time (which determines the adiabaticity of the process) could also be optimized for different values of J and noise parameters. However in this case we work with a single value for the sake of simplicity, which results in the fidelity saturating at a value that is smaller than 1.
To conclude, in this work we have presented a way to overcome two major bottlenecks in the quantum simulation of LGTs: one being the challenge of simulating fermionic matter and the other is the redundancy of the Hilbert space. We have shown how both of these problems can be tackled by solving the local constraints (Gauss' law) for the matter rather than for the gauge field.
To compare our method with the more conventional ones, we demonstrated it for the simplest case of Z 2 , in a one space dimension. In the conventional methods one uses the nonlocal Jordan-Wigner map; or solves the local constraints for  = 1). Each data point corresponds to a different value for the coherence time (assuming T 1 = T 2 ) and CZ gate duration, and the horizontal axis is the number of CZ gates that can be preformed in 10% of the coherence time (which is a measure of the quality of the quantum hardware). The dashed vertical line represents typical values for IBMQ machines. The inset in (a) shows the hardware requirements, defined as the number of CZ gates in 10% of the coherence time that is required for a ground-state fidelity of at least 90%; for different values of L. The blue squares are derived from the numerical data of (a), and the red triangles are linear extrapolation. the gauge field, which introduces non-locality; or both. Our method avoids both these techniques and thus extends easily to higher dimensions without non-locality, which we demonstrated for d = 2. Furthermore, the method can be applied to more complicated gauge groups, including non-Abelian ones (specifically U (N) and S U (N)), following the criteria worked out in Refs. [33,34], and used as a basis for quantum simulation, depending on the availability of suitable platforms (in terms of dimensionality and the Hilbert spaces required for the different gauge groups). In any case, the constraints and redundancy are eliminated, making the feasibility question technological and not conceptual.
Our scheme thus imposes fairly modest requirements on the simulator: no redundant components, no local constraints to maintain and no need to directly implement fermions. Moreover, we have shown how to implement it in one dimension using only simple single-and two-body operations, which we demonstrated experimentally on the IBMQ platform. While the platform is not suitable for implementing our method for the fully two-dimensional model, we demonstrated it on a quasi-two dimensional toy-model which is a minimal version of the theory for which traditional methods fail.
Because of these modest requirements we believe that our method could become an important and useful tool for quantum simulation of LGTs with fermionic matter. As quantum technology progresses, we expect that in the near future it will be applied to real unsolved problems rather than mere demonstrations.

ACKNOWLEDGMENTS
We acknowledge the use of IBM Quantum services for this work and to advanced services provided by the IBM Quantum Researchers Program. The views expressed are those of the authors, and do not reflect the official policy or position of IBM or the IBM Quantum team. This research was supported by the Israel Science Foundation (grant no. 523/20 and 2323/19).
G.P. and T.G. contributed equally to this work.

APPENDIX A: APPLYING THE MATTER REMOVAL TRANSFORMATION
Here, we give further details about the matter removal process, going from the hard-core bosonic setting to the one without matter at all. This should be read as a more detailed description of the process outlined in section III B.
It is straightforward to verify that the U (2) transformation (defined in section III B) gives rise to U (2) σ ± n U (2) † = P + n σ ∓ n + P − n σ ± n U (2) σ z n U (2) † = ε n S n σ z n U (2) X n U (2) † = σ x n X n σ x n+1 U (2) Z n U (2) † = Z n (68) which results in the transformed constraints of Eq. (31) simply by acting with U (2) on the hard-core bosonic version of Gauss' law (Eq. (25)). Acting with U (2) on the Hamiltonian, we find that the electric part is invariant, Considering the transformation of the mass Hamiltonian, we take a step back, and consider its effective form in the chosen sector. Using the relevant Gauss' law constraints (31) before applying U (2) , the mass Hamiltonian effectively takes the form As S n = Z n−1 Z n , this is invariant under the transformation U (2) , so we have the the effective expression for the transformed mass Hamiltonian (and "effective" here means "in the chosen sector"). Finally, consider the transformation of the interaction part. Here, we have terms like Z n−1 σ + n X n σ − n+1 (and its Hermitian conjugate, see Eq. (24)) that have to be transformed under U (2) , resulting in (using Eq. (68)) Noting that σ ± σ x = 1 2 (1 ± σ z ) is a projection operator to σ z = ±1 and that σ x σ ± = 1 2 (1 ∓ σ z ) is a projection operator to σ z = ∓1, we see that the transformed Hamiltonian is block-diagonal in the matter spins, with static σ z configurations. Having the constraints (31) in mind, we can restrict ourselves to physical states by ignoring all the terms but those that are projectors onto matter down-states (σ z = −1), and then neglect the matter spins altogether. Formally, Eq. (31) implies that where |out⟩ is a product state of all the matter spins in which they all point down, and ψ (2) is a state of the gauge fields. Then, we can define a Hamiltonian acting only on the gauge fields degrees of freedom, bỹ includingH (2) m = H (2) m ,H (2) E = H (2) E and H (2) GM = iJ n Z n−1 P + n X n P + n+1 + h.c.
= i J 4 n Z n−1 [X n + ε n (S n X n − X n S n+1 ) − S n X n S n+1 ] + h.c.
If we assume that J is real (for d = 1 this can be assumed without loss of generality), the expression will be simplified as only anti-Hermitian contributions in the sum would have to be considered. This eventually leads to the expression (36) in section III B.

APPENDIX B: TROTTER-ERROR ANALYSIS
We would like to approximate the time evolution under a Hamiltonian broken to three pieces, by the Trotterized sequence [25,26] e −iHt ≈ e −iϵH 1 e −iϵH 2 e −iϵH 3 N , where N = t/ϵ is a very large integer. The choice of N, as usual, is a compromise between the experimental capabilities and the Trotterization error that we allow. Suppose we allow some error δ. Then, we would like to have which can be bounded, in leading order, by [26] In our case, we have Clearly, the first and third commutator cancel each other, and thus the order choice presented in section III C, would be optimal. Upon computing the norms, we get that the error is of order implying that a reasonable N to use, given δ and the parameters h, j (independently of m) is N ∼ t 2 2δ J 2 + |Jh| L (83)

APPENDIX C: MEASURING NON-LOCAL OBSERVABLES
A valid question in the design of a quantum simulator is what are the observables one is interested in measuring, and how to measure them. In LGTs, the relevant quantities are the gauge invariant observables. In section III C we discuss the local ones, namely the number operator and the electric field operator, and explain how to measure them within our quantum simulation scheme. Here we focus on the the nonlocal ones, namely the mesonic string operators, for R ≥ 1 (R = 0 gives the local number operator). There are two questions to be asked; first, what are transformed opera-torsM (n, n + R), whose measurement in the simulated physical states ψ will correspond to measuring the original operators with respect to the original states |ψ⟩? In other words, what are theM (n, n + R) for which ψ M (n, n + R) ψ = ⟨ψ| M (n, n + R) |ψ⟩ .
The second question would be how to actually perform such measurements in our simulator.
The factor (−1) R(2n+R−1)/2 is irrelevant; and the mesonic string expectation value will be obtained from measuring the expectation value of the four termsM α (n, n + R), M 1 (n, n + R) = Z n−1 They are all products of X, Y, Z operators along the string, with spectrum ±1, and they can be measured, for example, by using an ancillary qubit interacting sequentially along the string with all the links, one after the other, and using the right controlled gates accumulating their X, Y or Z contribution to the product (this is not a new method -see, e.g., [76] for details).

APPENDIX D: TROTTERIZATION OF THE QUASI TWO-DIMENSIONAL HAMILTONIAN
We have to construct a Trotterization of Eq. (67) based on single-qubit rotations and two-qubit gates between neighbouring qubits on a three qubits chain. The non trivial terms are the three-qubit interactions Y 0 Z 1 Z 2 , Z 0 Z 1 Y 1 and the two-qubit terms Z 0 Y 1 , Y 1 Z 2 . Naively, we can implement each of these four terms using Eq. (44) and the analogous properties of the controlled-NOT (CX) gate, resulting in: Y 0 Z 1 Z 2 = CX 21 CZ 10 Y 0 CZ 10 CX 21 (89) Z 0 Z 1 Y 2 = CX 01 CZ 12 Y 2 CZ 12 CX 21 (90) where CX mn and CZ mn are controlled-NOT and controlled-Z operators between qubits m and n. This implementation requires 12 two-qubits gates per trotter step, so it is somewhat inefficient and we can do better. Instead of treating each of the four non-trivial terms separately, we note that we can implement two of them together, as, for example, one can verify that the transformation CZ 12 CX 21 CZ 10 takes Y 0 to Y 0 Z 1 Z 2 and Y 1 to Z 0 Y 1 . This means that we can implement the terms Y 0 Z 1 Z 2 and Z 0 Y 1 together, by transforming with CZ 12 CX 21 CZ 10 and rotating qubits 1 and 0 around the Y axis: e i J 2 (Y 0 Z 1 Z 2 +Z 0 Y 1 )t ≈ CZ 12 CX 21 CZ 10 U 0 Y (ϵ) U 1 Y (ϵ) CZ 10 CX 21 CZ 12 N , (93) where U n Y (ϵ) = exp (iϵ JY n /2), and a similar result holds for the other two terms: Z 0 Z 1 Y 2 + Y 1 Z 2 . These two steps together still use 12 two-qubit gates, but in this form the trotter step can be simplified further, since a combination of two entangling (two-qubit) gates on the same pair of qubits can be decomposed into a single entangling gate and a few single-qubit rotations. Using this kind of decompositions one can implement our trotter step with only 8 two-qubit gates, and 10 single qubit steps. The single-qubit depth can be larger if a general single-qubit rotation is not a native gate on the platform (as is the case for IBMQ devices), but on the other hand it is reasonable to assume that it can also be reduced with circuit optimization. We do not think it is useful to discuss this further here since two-qubit gates are the major source of error in current devices.