Digital quantum simulation of lattice gauge theories in three spatial dimensions

In the present work, we propose a scheme for digital formulation of lattice gauge theories with dynamical fermions in 3+1 dimensions. All interactions are obtained as a stroboscopic sequence of two-body interactions with an auxiliary system. This enables quantum simulations of lattice gauge theories where the magnetic four-body interactions arising in two and more spatial dimensions are obtained without the use of perturbation theory, thus resulting in stronger interactions compared with analogue approaches. The simulation scheme is applicable to lattice gauge theories with either compact or finite gauge groups. The required bounds on the digitization errors in lattice gauge theories, due to the sequential nature of the stroboscopic time evolution, are provided. Furthermore, an implementation of a lattice gauge theory with a non-abelian gauge group, the dihedral group $D_{3}$, is proposed employing the aforementioned simulation scheme using ultracold atoms in optical lattices.


Introduction
Gauge theories lie at the core of fundamental physics; the standard model of particle physics -describing electromagnetic, weak and strong interactions -is based on the principle of gauge invariance [1]. It requires introducing additional degrees of freedom, the gauge fields, to the matter fields: force carriers, mediating interactions between matter particles. If the coupling is small enough, perturbative expansions allow calculations up to arbitrary accuracy, as in QED (Quantum Electrodynamics). In some quantum field theories the coupling depends on the energy scale (running coupling) [2,3], and thus there are regimes where perturbation theory is not valid, e.g. QCD (Quantum Chromodynamics) at low energies. In such non-perturbative regimes only special methods can produce meaningful results. The most common approach so far has been lattice gauge theory [4,5]. The idea is to discretize space (or spacetime) to construct a framework in which numerical tools could be applied -with Monte Carlo methods being the most prominent ones [6]. In spite of their success (e.g. calculation of the low-energy hadronic spectrum of QCD [7]), there are limitations which are inherent to Monte Carlo simulations of lattice gauge theories. A major one is the sign problem, which prevents investigations in fermionic systems in finite chemical potential scenarios [8]. As a consequence, corresponding phases in quantum field theories still remain relatively unexplored, e.g. the quarkgluon plasma or the color-superconducting phase of QCD [9,10]. Another drawback of these simulations is that they take place in a Euclidean spacetime, thus making realtime dynamics inaccessible and preventing, for example, the study of non-equilibrium phenomena. One approach to overcome these obstacles is quantum simulation [11,12]. The idea is to build a highly controllable quantum system serving as a platform for simulations of another quantum system. In particular, quantum simulations of lattice gauge theories [13,14] have been proposed using various quantum devices, such as ultracold atoms in optical lattices [15][16][17], trapped ions [18,19] or superconducting qubits [20,21]. While the simulated models can be distinguished by features like the gauge group (abelian or non-abelian), the matter content (dynamical or static) or the dimension , there are also differences in the proposed simulation scheme. The first one -based on an idea of Feynman [47] -is to use a quantum computer (i.e. single and two qubit gates) to simulate the dynamics after Trotterization. The second one is an analogue approach: By appropriate engineering of the interactions, the Hamiltonian of the simulating system is exactly mapped to the desired one (which can be adiabatically changed), leading to an exact time evolution. The third one is a hybrid of both (e.g. [48]), where the time evolution is Trotterized but the different terms of the Hamiltonian are implemented using an analogue simulation, instead of quantum gates. It is important to note that the first simulation scheme will probably need quantum error correction, whereas the other two may not. Using the scheme suggested by Feynman, a trapped ion based quantum simulation of a lattice gauge theory was implemented in 2016 [49], allowing the observation of real-time dynamics in the Schwinger model, (1+1) dimensional QED. However, the simulation involved only four ions and it remains a big challenge to scale up such a system as it involves the construction of a quantum computer. In this work, we will focus on the third option. The main challenges of a quantum simulation of lattice gauge theories are threefold: First of all, to simulate dynamical matter, the simulating system must include fermionic degrees of freedom. Unlike in other quantum devices where fermionic statistics is imposed on spin degrees of freedom through Jordan-Wigner transformations, fermionic degrees of freedom occur naturally in ultracold atomic systems, as one can work directly with fermionic atomic species. This is beneficial in particular when dealing with two or more spatial dimensions. Second, gauge invariance, as the characteristic symmetry of lattice gauge theories, is not manifested naturally by the candidate quantum simulators. In analogue simulation schemes, where the degrees of freedom and the Hamiltonian of the investigated theory get exactly or approximately mapped onto the simulating system, local gauge invariance can be obtained either as a lowenergy effective symmetry [23,25,26] or by an exact mapping to an internal symmetry, like e.g. hyperfine angular momentum conservation [40,41]. Although the analogue approach works in one dimension (in particular as demonstrated by an ultracold atom experiment currently set up to study the Schwinger model [50]), it becomes problematic when considering the third requirement. The lattice gauge theory Hamiltonians in two or more spatial dimensions typically contain four-body interactions (the magnetic plaquette interactions). In the current analogue simulation schemes, this four-body term is realized only in fourth-order perturbation theory [40], thus leading to weak interactions and posing a major challenge on the way to higher dimensional quantum simulations of lattice gauge theories. This problem can be circumvented using the following concept: By introducing an auxiliary degree of freedom and entangling it with the physical degrees of freedom, the four-body interactions can be decomposed exactly as a sequence of simpler twobody interactions, resulting in stronger interactions compared to analogue simulation schemes. Because of the sequential nature of the entangling operations, during which all other interactions must be frozen, a stroboscopic time evolution is required. The time evolution is therefore decomposed into smaller pieces according to Trotter's formula: e −itH = lim N →∞ ( j e −itH j /N ) N [51]. This method has already been proposed in 2+1 dimensions to construct a digital scheme for lattice gauge theories with arbitrary gauge groups [31]. A concrete quantum simulation with ultracold atoms has been proposed for the groups Z 2 and Z 3 [31,32]. In this work we extend this proposal of an algorithm digitizing lattice gauge theories with arbitrary gauge groups to 3+1 dimensions. This is an important step towards the simulation of phenomena occurring in nature. To study the accuracy of the digital scheme, a thorough analysis of the digitization (Trotter) error is conducted. Another important goal is the simulation of gauge theories with non-abelian gauge groups. The second part of this work is therefore devoted to an ultracold atom based implementation of a lattice gauge theory with the simple non-abelian gauge group D 3 , following the general algorithm presented in the first part. This paper is organized as follows: First, a brief lattice gauge theory background will be provided, with an emphasis on the Hamiltonian formulation used later on for quantum simulation. In the second section the digital algorithm enabling quantum simulation of lattice gauge theories with dynamical fermions in three dimensions will be described. Afterwards, improved bounds on the digitization errors in lattice gauge theories will be given. In the last section, possible implementations based on ultracold atoms will be discussed, in particular the implementation of a lattice gauge theory with the dihedral gauge group D 3 , by exploiting its semidirect product structure.

Hamiltonian formulation of lattice gauge theories
Lattice gauge theories can be formulated in a Hamiltonian framework exhibiting a continuous time coordinate, as first proposed by Kogut and Susskind [52]. The lattice consists of d spatial dimensions, where the matter fields are placed on the vertices x ∈ Z d and the gauge fields reside on the links (x, k) (where k ∈ {1, .., d} denotes the direction in which the link points).
Since the matter particles are allowed to tunnel and thus their number is not conserved locally, the states on the vertices are described by elements of a fermionic Fock space. Assuming the gauge group G to be either compact or finite, we label its irreducible representations by j and represent the matter fields by spinors ψ † j m , where m denotes the components of j. Their behavior under group transformations, implemented by the unitary operator θ g , is (summing over repeated indices): where D j nm (g) is the irreducible unitary representation j of the group element g. We will work with staggered fermions [53], distributing the Lorentz components of the spinor over neighboring lattice sites such that occupied even sites will correspond to particles and vacant odd sites to anti-particles. The Dirac spinor is then regained in a continuum limit. The gauge transformationsθ g of staggered fermions are related to θ g by with e (o) denoting the even (odd) sublattice. We can define a state |D invariant under the above transformation (analogous to the Dirac sea in the continuum) where all odd sites are fully occupied and all even sites are vacant. The other physical ingredients, the gauge degrees of freedom, are described by a tensor product of local Hilbert spaces on the links. The elements of each single link Hilbert space can be expressed in the group element states {|g } g∈G . The group G can act on it in two ways, corresponding to left (L) and right (R) transformations: We define the group element operator U , a matrix of operators acting on the link Hilbert space: where for continuous groups dg is understood as the group (Haar) measure, whereas for discrete groups the integral reduces to a sum over the group elements. The hermitian conjugate of U in the Hilbert space and in matrix space are related by The group element operators obey the following rules under group transformations: (the j will be omitted in the following as only one fixed representation j is considered; generalization to more representations is straightforward). With these definitions at x + 1 x + 2 x + 3 hand we can define a local gauge transformation which acts on all degrees of freedom intersecting at a vertex. It depends on a group element g which itself can depend on the position (see Fig. 1 for illustration): where k is the unit vector in k-direction. A state |ψ is therefore said to be gaugeinvariant if Introducing the dual basis to the group element states, the representation basis {|jmn }, connected by the relation g|jmn = dim(j) |G| D mn (g) (with j labeling irreducible representations and m,n the components under left and right transformations), we can define a gauge-invariant "empty" state for the whole lattice, including matter and gauge fields: where |000 is a singlet state of the gauge fields in the representation basis, corresponding to the trivial representation. All other gauge invariant states can be obtained by acting with gauge invariant operators on this trivial state. A conventional lattice gauge theory Hamiltonian consists of four such types of terms: (i) The magnetic Hamiltonian One can obtain gauge invariant operators by taking products of U -operators along closed paths. The shortest such possible path is a plaquette, characterized by two directions k and l (k < l and l ∈ {2, .., d}). Adding over all pairs of k and l for every vertex x, one may construct: This term is called magnetic Hamiltonian as it corresponds to the magnetic energy in the continuum limit of the Yang-Mills cases.
(ii) The electric Hamiltonian The correspondence with the electric field becomes clear for the case of G = U (1) where -if we set f (j) = j 2 -the Hamiltonian is just a sum over the square of the electric field of all links. Similarly, for SU (2) f (j) = j(j + 1) corresponding to J 2 . The two terms above involve only gauge fields. They both add up to a generalized version of the Kogut-Susskind Hamiltonian for lattice gauge theories with compact gauge groups [54].
(iii) The fermionic mass Hamiltonian Introducing staggered fermions gives rise to the following staggered mass term: where the alternating sign comes from the Dirac sea picture: particles on even sites and anti-particles on odd sites.
(iv) The gauge-matter Hamiltonian The last term is a fermionic hopping term minimally coupled to the gauge fields in a gauge invariant way: The total Hamiltonian we want to simulate in the following chapters is the sum of all four pieces. The state defined in (9) is the non-interacting vacuum: the ground state of

Digital algorithm for the quantum simulation of lattice gauge theories in three dimensions
Interactions in typical quantum simulation platforms are usually two-body, e.g. atomic collisions in ultracold atomic setups or spin-spin interactions in trapped ion setups. Three-and four body processes are strongly suppressed on the relevant experimental timescales, making it much harder to map the Hamiltonian of the simulated model onto the simulating system, if the former includes interactions of more than two bodies. This is particularly relevant for lattice gauge theories since magnetic interactions are fourbody terms (see Sec. 2). For the purpose of quantum simulation of lattice gauge theories it is therefore desirable to design a scheme in which interactions involving three and more constituents can be rewritten as exact sequences of only two-body interactions. In this way, the energy scale associated to plaquette interactions is not limited by perturbative arguments (as in previous proposals) and the simulation can give access to a bigger region of the phase diagram. One approach to this problem is based on the idea of using an auxiliary degree of freedom that gets entangled with the physical degrees of freedom and mediates their interactions. In the following, we will briefly present an isometry which formalizes this idea (it is sometimes referred to as stator [55,56]). We anticipate that in this new framework the time evolution has to be realized stroboscopically due to the sequential nature of the entangling operations with the auxiliary system. Therefore, a digital algorithm based on Trotter's formula will be designed to simulate lattice gauge theories in three spatial dimensions, using only two-body interactions. This corresponds to the hybrid simulation scheme discussed in the introduction, where the time evolution is Trotterized but the individual parts of the Hamiltonian are still implemented by an analogue simulation. In the last section, bounds on the Trotter error will be provided.

Isometries
We consider two Hilbert spaces: H A representing the "physical" degrees of freedom, where the interaction is supposed to be implemented, and H B representing the auxiliary degrees of freedom (sometimes called control in the following This can be viewed as an entangling operation between the physical and the auxiliary degrees of freedom. If this entangling procedure is chosen in a certain way, operations on the physical Hilbert space can be implemented by acting only on the auxiliary state. Assume we want to realize a Hamiltonian H ∈ O(H A ) in the physical Hilbert space. For that, we need to create an isometry S and a hermitian operator H ∈ H B in the auxiliary Hilbert space in such a way that the following relation holds: An analogue relation for the time evolution follows directly, since H n S = SH n : Therefore, by creating such an isometry and acting with H on the control, we obtain the desired time evolution of the physical state |ψ A : The evolved physical state is still entangled with the auxiliary state which means that one can either perform another operation using the isometry S or disentangle both states. This would lead to a product state with the auxiliary state going back to its initial state:

The three-dimensional algorithm
In this section we discuss an algorithm to simulate the lattice gauge theory Hamiltonian in three spatial dimensions. We start from the lattice model described in Sec. 2. To create plaquette and gauge-matter interactions by means of isometries, we introduce an auxiliary degree of freedom in the middle of every second cube (either all even or odd ones) and assign to it a Hilbert space H isomorphic to the Hilbert spaces on the links (see Fig 2). Then, the lattice gauge theory Hamiltonian is split up into several parts which are implemented independently and sequentially: where we explicitly distinguish gauge-matter interactions taking place along different directions and in odd or even cubes, as well as plaquette interactions corresponding to the different plaquettes of a unit cube (therefore we get a sum of six terms in both cases). The desired time evolution e −itH LGT is then approximated by a Trotterized time evolution consisting of N steps: e −itH LGT ∼ ( j e −itH j /N ) N , where H j is any of the terms appearing in (20). While electric and mass terms can be treated easily using only the physical degrees of freedom, the plaquette and gauge-matter terms are further decomposed as a suitably chosen sequence of simpler interactions mediated by the auxiliary systems. This sequence will then be executed in parallel for all cubes where auxiliary degrees of freedom are located. However, since for the gauge-matter interactions the individual parts of this sequence do not commute for adjacent links, we have to place the auxiliary d.o.f. in every second cube to avoid undesired interactions. The exact decompositions will be given in the next sections.

Plaquette interactions
Since we put auxiliary atoms in every second cube, we can not realize all plaquette interactions at once and we split them up in the following way: It is important to mention that the six magnetic terms commute, therefore e −iτ H B = j e −iτ H B,je e −iτ H B,jo and this splitting does not affect the error of the Trotter approximation (20). To implement each term we will use the isometry where the first Hilbert space belongs to the gauge field, residing on link i, and the second one to the aforementioned auxiliary degree of freedom in the center of the cube. It fulfills the relation allowing to realize operations on the link i through the auxiliary degree of freedom. The isometry S i can be created by the unitary acting on the initial state | in = |ẽ . We repeat similar entangling operations U i (or U † i ) for the three other links of the plaquette (e.g. the links 1,2,3,4 of cube x, see Fig. 3) and obtain a plaquette isometry of the form The crucial part is that it fulfills the relation i.e. acting locally with on the control of cube x enables us to realize the magnetic time evolution for this plaquette. The required sequence acting on the plaquette state |ψ 1234 , the tensor product of the four link states, and the auxiliary state | in is The other two plaquette terms associated to cube x can be created in the same manner but with different isometries. Using the abbreviations for the gauge field operators defined according to Fig. 3, we need to replace the isometry S 1234 (x) by Applying the sequence from (28) gives then rise to the time evolution of the physical state under H B,2 (x), or H B,3 (x).
We can now formulate an algorithm to implement the whole plaquette interactions. We start with the controls placed in the center of every even cube and do the following three steps: (i) Create the isometry: Let all the controls interact with all the gauge fields on links of type 4 and create the unitary . Repeat similar processes with links 3, 2 and 1 to obtain the unitaries x even U 1 (x). In total, we get: x even (ii) Act on the controls with the Hamiltonian (iii) In the last step, undo the isometry by creating the inverse of the first step, i.e.
x even The above procedure is applied to a state |ψ | in . Thanks to relation (28) we obtain: x even We repeat the procedure with the two isometries S 5671 and S 5894 . In this way we create W B,2e = e −iH B,2e τ and W B,3e = e −iH B,3e τ . The same steps are then repeated with the auxiliary degrees of freedom moved to the center of the odd cubes so that we can implement W B,1o , W B,2o , W B,3o . Since all pieces of the magnetic Hamiltonian commute, this sequence gives us exactly the magnetic time evolution:

Gauge-Matter interactions
After expressing the four-body plaquette interactions as a sequence of two-body interactions, we want to obtain the gauge-matter interactions in a similar way. We need again to split up the relevant Hamiltonian terms into parts suitable for implementation: An important ingredient for rewriting these interactions as two-body terms is the following unitary operation, entangling the fermion at vertex x and the gauge field on link (x, k): where Z mn = −i(log mat (U )) mn , and the logarithm is taken only in matrix space (welldefined since the matrix elements commute). Its meaning becomes more apparent if we assume the gauge group G to be compact; then, we obtain an interaction of the "vector potential" operatorφ a with the fermionic charge Q a . It can therefore be interpreted as a fermionic transformation whose parameter is an operator acting on the gauge field. The idea is now to use this transformation to map a pure fermionic tunneling term into the desired gauge-matter interactions, as Thus, defining the fermionic tunneling Hamiltonian as allows writing the Hamiltonian H GM as: Since every fermion is connected to six links in three dimensions we have to split up the process in six steps as described in the beginning. We start by realizing H GM,1e , i.e. H GM (x, 1) for all even links (see Fig. 4). We apply the following sequence: (i) Let the gauge degrees of freedom interact with the fermions at the beginning of the link to obtain the unitary: x even U † W (x, 1) .
(iii) Let the link degrees interact again with the fermions to generate: x even

This gives us in total
x even By applying a similar sequence for the other links of the cube, we can create W GM,2e , W GM,3e , W GM,1o , W GM,2o , W GM,3o .
Using isometries, there is an alternative way of realizing the gauge-matter interactions. It requires more steps but on the other hand does not require interactions between the physical degrees of freedom as all of them are mediated by the auxiliary degrees of freedom. The sequence goes as follows: (i) Let the controls -initially placed in all even cubes in the state | in = |ẽ -interact with the gauge links U 1 according to (24) to create the isometry S 1 : x even (ii) Let the control interact with the fermion at vertex x to realize the interaction x even U † W (x, 1) which is the same interaction as U † W (x, 1) but between the control and the fermion ψ m (x). Due to the properties of the isometry S 1 the interaction between the control and the fermion will translate into an interaction between the fermion and the link.
(iii) Afterwards, allow for pure tunneling between the fermions which gives rise to x even e −iHt(x,1)τ .
(iv) Following (35), apply U W (x, 1) for all even cubes which is again realized by an interaction between the control and the fermion ψ m (x): x even U W (x, 1) .
(v) Finally, we have to undo the isometry between the control and the gauge field: x even The resulting sequence -applied to some physical state |ψ and the auxiliary state | in -is: x even We repeat a similar procedure for all other links in the cube which gives us W GM,2e , W GM,3e , W GM,1o , W GM,2o , W GM,3o .

Other parts of the Hamiltonian
The electric part W E = e −iH E τ and the matter part W M = e −iH M τ are local terms of our Hamiltonian and thus one can implement them by acting locally on the physical degrees of freedom. We can now write down the whole sequence for a time step τ (combining commuting magnetic terms to W B ): It's important to notice that all time evolutions in the above sequence are individually gauge-invariant. Therefore, errors coming from the digitization do not break gaugeinvariance.

Error bounds for Trotterized time evolutions in lattice gauge theory
Although the approximated Trotter evolution has the correct gauge symmetry, it is still important to analyze how much it deviates from the desired exact time evolution.
In this section we derive bounds for the Trotter error, according to the digitization scheme presented in the previous section. We focus on the standard Trotter formula (the first order formula) and the second order formula which gives a better approximation without major changes in the implementation. We do not consider higher order formulas, because they would require more experimental effort in the sense that the tunability of the experimental parameters would have to be much more flexible and the number of operations required for a single time step would increase exponentially with the order of the approximation [57]. The first order formula [51] is of the form For which, using the operator norm, the difference to the physical time evolution U(t) = e −itH can be bounded by [58][59][60]: To get a better scaling with the number of time steps we can apply the Trotterization sequence in reverse order after the usual Trotterized time evolution (second order formula) [61]: From an implementation point of view this decomposition can be realized straightforwardly once we know how to obtain the sequence for the first order. Following the proof in [62] adapted to unitary operators, an upper bound for the trotter error can be derived: Compared to the first order formula, the second order formula has an error which decreases faster with the number of time steps N at the cost of a longer sequence. The experimental difficulty, however, is the same for both decompositions. We can now specify these bounds for lattice gauge theories. This is an important task since an implementation of this digital scheme will have to balance experimental errors, which can break gauge-invariance and increase with the number of steps in the sequence, and errors caused by the digitization, which have the opposite behavior. Therefore, a precise bound of the Trotter error helps in finding the optimal number of steps, so that experimental errors do not accumulate unnecessarily and the chance of breaking gauge invariance is reduced as much as possible.
Since the different parts of the Hamiltonian can not be implemented simultaneously, they are split up in the digitized simulation scheme. Hence, for the computation of the trotter error we divide the Hamiltonian into these individual pieces, according to the Trotterized time evolution given in (38). Generalizing to d dimensions: 3.3.1. First order formula By inspection of (40) we see that for an upper bound on the digitization error of the standard trotter formula, the commutators among all different parts of the Hamiltonian in (43) have to be evaluated, as well as their norms. Since the derivations are very lengthy we will refer the interested reader to the Appendix. We provide here the final result: where d is the number of spatial dimensions, d U the dimension of the representation of the group element operator U and N links the number of links in the lattice. One might think that operator norms involving H E are unbounded but, since we either work with finite groups (whose number of irreducible representations is finite) or appropriate truncations of infinite gauge field Hilbert spaces, the expression max j |f (j)| is finite, so that we always obtain sensible error bounds.

Second order formula
To bound the error of the second order formula we need to calculate nested commutators according to (42). Details on their calculation can be found in the Appendix. We provide here the final result: If we assume a cubic lattice with L lattice sites per side we can express the number of links as: The upper bound shows that N should scale as N ∼ L d/2 t 3 which is somewhat bad since it considers a very general setting. If we restrict ourselves to the observation of intensive quantities we expect this scaling to be much better. However, there are observables in lattice gauge theories, e.g. Wilson loops, which do not fulfill this requirement and thus need to be bounded by more general estimates like the ones given above.

Implementation of digital lattice gauge theories with ultracold atoms
With this general scheme for the digital construction of three-dimensional lattice gauge theories at hand, we can turn to the implementation of some concrete examples with ultracold atoms. Typical gauge groups of interest are compact (e.g. U (1)), for which the link Hilbert spaces are infinite. A truncation of this Hilbert space is therefore required to make the quantum simulation feasible. Previous proposals have performed this truncation in the representation basis [13,14]. This procedure, however, spoils unitarity of the group element operators U and prevents the use of isometries (see 3.1). Thus, the Hilbert space of the gauge field should be truncated using group element states instead. A truncation of U (1) in this sense is given by the finite groups Z N which converge to U (1) in the N → ∞ limit. The digital quantum simulation of Z N gauge theories has been studied in [31,32]. We summarize below their main features, and then we build on these to tackle the simulation of simple non-abelian gauge models with dihedral symmetry given by the group D N .

Implementation of lattice gauge theories with gauge group Z N
Lattice gauge theories with a finite abelian gauge group play an important role as they approximate compact quantum electrodynamics. Since the Hilbert space of the gauge field is reduced to dimension N if the gauge group Z N is considered, ultracold atoms can be used to represent these gauge degrees of freedom. These N states are labeled by |m and we define unitary operators P and Q on them: Since the group is abelian, its representations are one dimensional and we need to consider a single fermionic species, ψ † , on the vertices. We can now define the Hamiltonian of Z N lattice gauge theory with fermionic matter: Possible implementations for Z 2 [32] and Z 3 [31] with isometries have been discussed in two space dimensions. These proposals can be readily generalized to three dimensions following the scheme presented in the previous section. The matter content is represented by a fermionic atomic species whereas the gauge fields can be represented by a second atomic species with the appropriate ground state manifold, e.g. F = 1/2 for Z 2 or F = 1 for Z 3 . Furthermore, auxiliary atoms must be trapped in the center of each second cube. These species are confined to the desired lattice geometry by suitable optical lattices and their interactions are realized by ultracold atomic scattering.
Since the type of interactions appearing in two and three dimensions are the same, the implementation in three dimension follows closely the steps explained in [31,32] and the reader should refer to the original references for more details.
Here we just report the bounds on the Trotter error that can be computed following the discussion in Sec. 3.3. In three dimensions and for the gauge group Z N , we obtain the first order formula (see (44)) : and the second order formula (see (45)): Note that these formulas give a more accurate bound with respect to the original analysis in [31,32].

Implementation of lattice gauge theories with a dihedral gauge group
We now turn our attention to the implementation of simple non-abelian lattice gauge theories, with symmetry given by the dihedral group D N (with N odd and N ≥ 3 which converges in the large-N limit to O (2)). This symmetry group can be characterized by a set of rotations R in a two-dimensional plane and reflections S along a certain axis: The above notation already suggests that D N can be decomposed into a semidirect product of the abelian groups Z N and Z 2 corresponding to rotations and reflections: It is thus useful to write the states of the gauge field Hilbert space as states living in the tensor product of an N -dimensional Hilbert space and a twodimensional one, |p, m = |p ⊗ |m ∈ H N ⊗ H 2 . In the implementation, such a product Hilbert space can be realized by using two atoms with the appropriate hyperfine structure. If we choose to work with the smallest faithful irreducible representation of the group, we need two different fermionic components for the matter, denoted by ψ 1 and ψ 2 , due to the non-abelian nature. Accordingly, the gauge field operators U on the links become matrices of operators U = e i 2π Np σz σm x (p acts on H N andm on H 2 ; σ x and σ z act in matrix space). This allows us to write down the Hamiltonians where f r and f l satisfy the condition f l = f −l ∀l. D N lattice gauge theories do not have a meaningful large-N limit (like Z N with compact QED) as O(2) is not a "conventional" lattice gauge theory and does not have a continuum limit. Thus, in principle the coefficients in (52) can be chosen arbitrarily. However, it is convenient to identify the second term of (52) with the electric energy of a Z N lattice gauge theory (see above), and fix the coefficients accordingly.

Simulating system
Our implementation scheme is in principle applicable to all dihedral groups but we focus here on the simplest case D 3 (isomorphic to the group of permutations S 3 ). We first discuss the system we will use as a platform to perform the quantum simulation. For the simulation of the matter fields it is crucial to use fermionic atoms to obtain the correct commutation relations. A natural, minimal choice for the two fermionic d.o.f. ψ 1 and ψ 2 is to use the two internal levels of an atom with a F = 1/2 hyperfine ground state. For example ψ 1 and ψ 2 can be associated with the F = 1/2 multiplet in the following way: These atoms must be trapped by a superlattice that allows to modulate the depth of the minima (to account for the staggering) and the tunneling rate between nearest neighbors (to switch tunneling on and off in the different steps of the Trotter sequence).
To simulate the gauge field and auxiliary Hilbert spaces, we will exploit the product structure as mentioned above: H aux H link H 3 ⊗ H 2 . One convenient choice is to use two atomic species: a bosonic one with an F 3 = 1 hyperfine multiplet (the index 3 will label the three-level system) and a fermionic one with an F 2 = 1/2 multiplet (the index 2 will label the two-level system). In total, we need four different atomic species: two atoms trapped at the middle of each link, and two extra atoms (that must be addressed independently of the previous two) in the middle of each second cube. For the links, we identify: Every state of the Hilbert space on the link can be obtained as a tensor product of the two multiplets, e.g. |p = 1, m = 1 = |F 3 = 1, m F = 1 ⊗ |F 2 = 1 2 , m F = − 1 2 . The corresponding creation operators on some link (x, k) are described by a † m F (x, k) with m F = −1, 0, 1 for the three-level system and c † m F (x, k) with m F = −1/2, 1/2 for the two-level system. It is useful to introduce unitary operators P 3 , Q 3 and P 2 , Q 2 acting respectively on the three-level and two-level atoms. They are defined as: The operators P 3 , Q 3 fulfill the Z 3 algebra whereas the operators P 2 , Q 2 fulfill the Z 2 algebra. The Hilbert space of the auxiliary atoms has the same structure, and we label its states/operators with a tilde to distinguish them from the corresponding link quantities, i.e. we have states |p and |m and operatorsã † m F (x) (with m F = −1, 0, 1) andc † m F (x) (with m F = −1/2, 1/2). The link and auxiliary atoms must be trapped in the desired positions by arranging suitable optical potentials. The individual minima must contain exactly one atom and must be deep and well separated so that the dynamics is frozen (no tunneling, no interactions between nearest neighbors). When requested, the lattices must undergo a rigid translation so that specific pairs of atoms can overlap and interact via twobody scattering. The resulting setup -for convenience projected to two dimensions -is depicted in Fig. 5. All interactions between the constituents of the simulating system from above are in the form of two-body scattering. As will become clear in the following, we need to impose specific constraints on the scattering. First we want interactions that are diagonal in m F and do not change the internal level of the atoms. This can be achieved by lifting the degeneracy of the hyperfine multiplets such that transitions changing m F will cost energy. A possible way to do this is by introducing a uniform magnetic field which adds the following correction to the Hamiltonian (Zeeman shift): Figure 5. The simulating system consists of one atomic species on the vertices representing the matter (red) and two for both the gauge fields (blue) on the links and the controls (green) located at the center of every second cube (projected into two dimensions for better visualization). The simulated degrees of freedom are encoded in the hyperfine structure of the atoms, i.e. either an F = 1 or an F = 1/2 multiplet. The alternating occupation of vertices with fermionic atoms shall illustrate the staggered fermion picture, in which this configuration corresponds to the non-interacting vacuum (see Dirac sea in the continuum). The empty green circles indicate the need to move the auxiliary atoms between even and odd cubes.
where µ B is the Bohr magneton and g F the hyperfine Lande factor. The energy splitting has to be different for different atomic species to avoid resonant exchanges, therefore we need to choose species with different Lande factors. Another possible approach to realize the different energy splittings is to address each species individually, for example exploiting the AC Stark effect. Second, at some point we need to modulate the interaction strengths depending on the internal level of the atoms. This can be achieved for example by spatially separating the different m F levels via a magnetic field gradient. The different m F levels will experience forces pointing in different directions and reach different equilibrium positions within the same potential well. By properly choosing the Lande factors of the atomic species and tuning the magnetic field gradient one can then tailor the overlap of the atomic Wannier wave functions (and hence their interaction strength) in an m F -dependent way. Below we discuss several details of the implementation scheme.

Initial configuration and background Hamiltonian
Before starting the simulation we should define the initial configuration of our simulating system. It is reached if all optical potentials are sufficiently deep and separated such that no tunneling occurs and all atomic wave functions do not overlap. All minima of the auxiliary lattice are loaded with one atom in the group element state corresponding to the identity group element, i.e | in = |ẽ = |0,0 . This means we have to prepare the state | F 3 = 1;m F = 0 for the three-level system and | F 2 = 1/2;m F = 1/2 for the two level system. The preparation of the atoms representing gauge and matter fields depends on the initial physical state we want to simulate. All atoms must occupy the motional ground state with energy E 0 (different for different atomic species). As mentioned in the previous section, we also introduce a uniform magnetic field (or an AC Stark effect) to lift the degeneracy of the ground state manifolds and induce energy splittings ∆E (again different for different species) between the different m F components. We can define the non-interacting Hamiltonian H 0 which will be present throughout the whole implementation: All parts of the digital simulations are added on top of H 0 . To recover the desired Hamiltonian H of our D 3 lattice gauge theory, we move to an interaction picture , i.e. we will work in a rotating frame with respect to H 0 and make use of the rotating wave approximation.

The mass Hamiltonian The mass Hamiltonian in three dimensions takes the form
with ψ † (x)ψ(x) = ψ † 1 (x)ψ 1 (x) + ψ † 2 (x)ψ 2 (x). Thus, the corresponding time evolution W M = e −iH M τ for a time step τ can be implemented by smoothly modulating the superlattice trapping the fermions so that the energy of the even minima is increased by an amount M even . This results in the Hamiltonian If we act with this Hamiltonian for time Meven M τ , we obtain the desired unitary evolution W M , up to an irrelevant global phase.

Creating the isometry
The creation of plaquette interactions and gauge-matter interactions involves constructing the isometry S i (see Sec. 3), entangling auxiliary atoms with the atoms on link i. If we want to create it from the auxiliary state corresponding to the neutral element | in = |ẽ , we have to apply U i = dg |g i g| i ⊗Θ L † g . Specifying this equation to the gauge group D 3 , we obtain the following interaction between the d.o.f. on link i and the ones of the control: wherep = p p |p i p| i andm = m m |m i m| i . As defined previously, Q 2 and Q 3 are the raising operators of the auxiliary atoms, i.e Qm 2,i and Qp 3,i raise them F -values of the auxiliary atoms according to the m F -values of the atoms on link i. By choosing |0,0 as the initial state of the auxiliary atoms, the creation of the isometry reduces to an interaction between the three-level atom on the link and the auxiliary three-level atom in parallel with an interaction between the two-level atom on the link and the auxiliary two-level atom. These are the same interactions required for creating the isometry of a Z 3 lattice gauge theory [31], respectively a Z 2 lattice gauge theory [32]. We can therefore directly adopt the procedure from [31,32]. The idea is to write Qp 3,i and Qm 2,i as an interaction between the z-components of the hyperfine angular momentum operators F z,3 and F z,3 , respectively F z,2 and F z,2 : where V † 3 is a local change of basis from the P 3 -basis {|p } to its conjugate Q 3 -basis and: where V † 2 is mapping from the P 2 -basis {|m } into the conjugate Q 2 -basis. The basis transformations V 3 and V 2 are local operations on the auxiliary atoms that can be implemented with optical/RF fields. The interactions between the z-components of the hyperfine angular momentum operator can be realized by introducing an energy splitting between the different m F -levels such that the two-body scattering term will contain only m F preserving terms. The sequence to obtain U i is therefore: To undo the isometry it is necessary to create the conjugate of these interactions which can be done by flipping locally them F = 1 andm F = −1 levels, thus mapping F z,3 into − F z,3 . In the same way, the conjugate of the two-level system is created.

4.2.5.
Plaquette interactions Knowing how to construct the isometry, the implementation of the plaquette interactions is straightforward. Since we have to split them in six different parts (see Sec. 3), we start with H B,1e , the type 1 plaquettes of the even cubes, where the auxiliary atoms are placed in the standard configuration. We follow the three steps of the algorithm given in Sec. 3.2.1: (i) We create the plaquette isometry out of the isometries S i which is realized for a link i by moving the lattice of the auxiliary atoms to the respective link and tailoring the interactions as described above (neglecting the basis transformations V for the moment). This can be done in parallel for the whole lattice: The desired plaquette isometry is obtained by applying this procedure to all four links and including overall basis transformations V 3,all and V 2,all : x even This operation, acting on the tensor product of |0,0 and any state of the links, gives rise to the proper entangled state which maps plaquette interactions to local operations on the control.
(ii) The next step is a local operation on the auxiliary Hilbert space. We need to implement e −i H B τ with H B being the control Hamiltonian H B = λ B Tr( U + U † ). This requires an interaction between the two-level and the three-level system: wherem = m m |m m|. We can rewrite P 3 + P † 3 = −I + 3 |0 0 | ≡ −I + 3N 0 with N 0 ≡ã † 0ã 0 . Defining a number operator for the | F 2 = 1/2; m F = 1/2 state of the two-level system as N 1/2 ≡c † 1/2c 1/2 we can write down the interaction e −i H B τ : The first exponential is a local term of the two-level system which can be implemented by means of optical/RF fields. The second term requires scattering between the two auxiliary atoms. The corresponding Hamiltonian density in second quantized form is [63]: where Φ † α denotes the creation operator of the atomic Wannier wave function corresponding to the internal state α and µ the reduced mass of the two atomic species. The projection operators onto the different scattering channels are expressed by polynomials of F 1 · F 2 , the coefficients g k are therefore functions of the scattering lengths. To obtain the time evolution due to this interaction we have to integrate the Hamiltonian density over space and time. Since eq. (68) involves only specific levels, we need to turn on the magnetic field gradient and split the different m F components such that only them F = 0-component and them F = 1/2component overlap during the collision. Moreover, changes in the internal level of the two atoms during the collision are suppressed by the Zeeman splitting. With these assumptions, the time evolution is described by the following unitary with g 0 = 1 6 (3a 1/2 + 4a 3/2 ) (a 1/2 , a 3/2 are the scattering lengths for the scattering channels with F tot = 1/2 and F tot = 3/2) and α the time-integral of the overlap of the two wave-functions during the collision. By carefully tuning the interaction time we can set α = 6λ B τ g 0 and finally obtain: which is up to local operations the desired unitary V B . This interaction will be implemented in parallel for all cubes where auxiliary atoms are placed, i.e. in this case for the even cubes. Hence, the overall interaction of this step is e −i x even H B (x)τ . When the magnetic field gradient is on, different levels of the hyperfine multiplet will acquire an extra energy splitting with respect to the background Hamiltonian (57). This induces extra phases that need to be cancelled somehow. For example, after the collision has been completed, we can invert the slope of the gradient and accumulate phases in the opposite direction until the net effect is zero (this trick has to be applied for all scattering events of this kind).
(iii) In the third and last step we have to undo the isometry. This can be done by taking the hermitian conjugate of the first step, i.e. the sequence: According to (29) these three steps give us W B,1e . If we repeat now the same procedure but with the links corresponding to the second and third plaquette term, we obtain W B,2e and W B,3e . To realize the odd cubes time evolution, we move the auxiliary atoms to the centers of the odd cubes and repeat all of the above. This results in the time evolutions W B,1o ,W B,2o and W B,3o . Afterwards, the auxiliary atoms are brought back to the centers of the even cubes.

Gauge-matter interactions
For the Gauge-matter interactions on a link (x, k) we have to implement the Hamiltonian with U p ≡ e i 2π 3 σzp and U m ≡ σm x . We can use the product structure of U to implement the gauge-matter part via two-body interactions. We follow the procedure given in Sec. 3.2.2 and define the unitaries U W , one corresponding to U p : and another one corresponding to U m : (74) With these definitions at hand we can get the following relation by applying twice the Baker-Campbell-Hausdorff formula: The gauge-matter Hamiltonian can then be written as with the tunneling The crucial thing to note here is that all the terms involve only two-body interactions which allows an implementation with the proposed ultracold atomic setup. We can not implement all gauge-matter interactions at once as the fermions on the vertices are only allowed to interact with one link at a time. Focusing on the links in the 1-direction for the even cubes, we describe how to realize the time evolution e −i x even H GM (x,1)τ . Since we want to keep the lattice of the matter and link degrees of freedom fixed, these interactions will be mediated by the control atoms according to the algorithm presented in 3.2.2.
(i) We first build the isometry S 1 between auxiliary atoms located at the center of even cubes x and the corresponding atoms on link (x, 1), x even U 1 (x). This interaction can be implemented exactly in the same way as already done for the plaquette term (see (64)). Due to the relation in (23) the gauge-matter interactions will then translate into an interaction of exactly the same form but between the auxiliary atoms and the fermions.
(ii) Afterwards, the two terms U † W,p and U † W,m have to be implemented by two-body scattering processes but between the fermions and the auxiliary atoms due to the isometry, therefore denoted as U W,p and U W,m . Starting with U † W,p , we first write it in terms of the angular momentum operator respectively the second quantized operators ψ 1 and ψ 2 for the fermions: Now we have to tailor the atomic collision between the F 3 = 1 and the F = 1/2 multiplet accordingly. The magnetic field again lifts the degeneracy of the hyperfine levels and thereby prevents any transitions changing the m F -values. The interaction Hamiltonian contains two possible scattering channels and gives rise to the following time evolution: with g 0 = 1 6 (3a 1/2 + 4a 3/2 ), g 1 = 2 3 (a 3/2 − a 1/2 ) (a 1/2 , a 3/2 are the scattering lengths for the scattering channels with F tot = 1/2 and F tot = 3/2) and β the time-integral of the wave-function overlap. If we tune overlap and interaction time such that β = 2π 3g 1 we obtain The second exponential is the desired interaction U † W,p whereas the first term is a fermion-dependent phase, denoted from now on as where θ = 2πg 0 3g 1 and ψ † ψ = ψ † 1 ψ 1 + ψ † 2 ψ 2 . A discussion of these phases will be done later on. Before, the implementation of U † W,m is explained. It has the form: with N −1/2 ≡c † −1/2c −1/2 and V H,f er = 1 √ 2 (σ x,f er + σ z,f er ) a Hadamard transform on the fermions which can be implemented by means of optical/RF fields. The remaining two-body interaction is realized as scattering between the F = 1/2 states of the control atoms and the fermions. It can be described by the following unitary: (for the explicit form of g k see [32]). We switch on a magnetic field gradient designed in a way that only the m F = −1/2 -components of the auxiliary atom and the fermion overlap. Moreover, the interaction time should be tuned such that γ = π g 0 +g 1 which gives rise to: Since the implementation of U † W,p and U † W,m is done in parallel for all even cubes we get the sequence (iii) In the next step we implement the tunneling in the 1-direction for even cubes which can be achieved by modulating the superlattice and decreasing the potential barriers on the desired links. We get x even e −iHt(x,1)τ (85) (iv) After the tunneling we need to realize the conjugate of U † W,p and U † W,m , i.e. U W,p and U W,m . One way of creating U W,p is by doing a spin flipping operation V F,3 for the three-level system of the control which results in: For the creation of U W,m we simply observe that U † W,m is hermitian. The sequence for step 4 is (v) In the last step we need to undo the isometry, which is done by the conjugate of the first step, x even We summarize by writing down the whole sequence acting on the initial auxiliary state | in = |0,0 : x even x even (88) We finally get the desired gauge-matter interactions up to the fermionic phases V W (θ). However, if we consider the whole lattice (on which the number of fermions is globally conserved) it can be shown that the phases correspond to a static vector potential of zero magnetic field and are therefore unphysical, as carried out in the procedure given in [31]. If we repeat the whole sequence (88) for the other links we obtain the gaugematter interactions W GM,2e , W GM,3e , W GM,1o , W GM,2o and W GM,3o .

Electric Hamiltonian
The electric Hamiltonian for the gauge group D 3 acts on the gauge fields residing on the links. If we choose its second part -which corresponds to pure rotations only -in accordance with the electric energy of Z 3 we obtain, using the notation of previous sections: If we also express the interactions of the first part in terms of operators acting on the link atoms, we end up with: The first Hilbert space represents the three-level system, the second one the two-level system. The coefficient f l is the overall coefficient for the electric part corresponding to pure rotations, equivalently to Z 3 . We have to implement the time evolution: The first and the third exponential are local terms of the atoms and can be addressed by external fields. The second term is implemented by two-body scattering similar to the one for the plaquette interactions. Therefore, we need to bring the two atoms together, which should be simple to implement since both of them are trapped near the middle of the link. Following the steps for the plaquette interactions, we obtain: Tuning overlap and interaction time such that δ = λ E frτ g 0 and combining it with the local operation V 2 = e i λ E fr τ 2 N 1/2 , gives us: If we then perform a Hadamard transform V H,2 on the two-level system, we get the desired interaction: which gives us the electric Hamiltonian up to local operations. We have implemented all interactions using local operations on the atoms and tailoring the appropriate two-body scattering terms. If we use the sequence to evolve the system for a time τ = T /N and we repeat the same sequence N times, we get a Trotter approximation of the desired time-evolution e −iH LGT T . The accuracy of this approximation is discussed below.

Errors
The errors affecting the precision of the simulation are twofold. On the one hand, we have Trotter errors coming from the digitization which can be estimated by specifying the general error bounds given in Sec. 3 to the case of three dimension and gauge group D 3 . We obtain for the fist order formula (see (44)): and the second order formula (see (45)): We stress again that the digitization error doesn't break gauge invariance, because all steps of the sequence individually respect the right symmetry. Therefore, the Trotter expansion can only give rise to quantitative deviations, but not to qualitative changes. On the other hand, there will be experimental errors in the implementation. Unlike errors caused by the Trotterization, they may break the gauge symmetry and accumulate step by step. We briefly want to look at the scaling of these errors. We consider a small perturbation h j to one of the Hamiltonians H j which is realized during the implementation of the Trotter sequence (38). The difference of the time evolution e −iτ (H j +h j ) to the desired one e −iτ H j can be bounded to first order in the operator norm by h j τ . To get the total experimental error caused by the gates corresponding to H j , we need to look at the whole Trotterized time evolution (we focus here on the second order formula (41), i.e. the gate is repeated 2N times). We have to distinguish four cases: On the one hand, whether the experimental error is statistical or systematic and, on the other hand, whether the implemented gate depends on the simulated time (e.g. electric Hamiltonian, fermionic tunneling, etc.) or not (e.g. entangling operations). The advantage of a statistical error is that we can apply the central limit theorem and obtain a scaling of √ N with the number of Trotter steps compared to a linear scaling in the case of a systematic error. In the same vein, a gate depending on the simulated time t is advantageous since the time step τ = t 2N in each trotter sequence scales as 1 N , whereas for gates not depending on t, the error scales with some fixed amount of time t exp,j specific to the gate. The bounds for these four types of experimental errors are summarized in Table 1. We see that operations that do not depend on the simulated time t are the ones most prone to errors. During their implementation a lot of care should be taken, in particular to avoid systematic errors. When estimating the error of the whole implementation sequence, one should keep in mind that errors of different gates are generally independent and thus do not add up linearly. However, the total experimental error will still increase with N , so that the number of Trotter steps has to be chosen in a way to balance digitization and implementation errors.
Typical sources of errors in ultracold atom experiments are as follows: The first one is decoherence, e.g. caused by spontaneous scattering of lattice photons with the atoms, atomic collisions with the background gas, field (laser or magnetic) fluctuations, etc. This is relatively well under control nowadays, where coherence times t coh of the order of minutes have already been achieved [15,64,65], thus requiring the total simulation time t sim to fulfill t sim << t coh . Secondly, one needs to ensure that the atoms remain in the lowest Bloch band throughout the whole implementation. Hence, it is of crucial importance to shape the lattice and move the atoms in an adiabatic way. This is particularly important in our simulation scheme, where the auxiliary atoms have to be moved around or when the matter lattice has to be deformed to allow tunneling. This means that the corresponding timescale t mov should be bigger than the inverse of the frequency ω associated to the energy difference between lowest and first excited Bloch band (t mov >> 1/ω), while at the same time the obvious constraint t mov << t coh has to be fulfilled. However, such techniques have also become well-controlled [15,66]. Errors more specific to this proposal are connected with the tailoring of the twobody scattering. This requires a high degree of control over the overlap of the atomic wave functions and accurate timing of interaction during these collisions. This is also dependent on the ability to design and manipulate the magnetic field gradient in a precise manner.

Summary
In this work, two main results were discussed. First, a digital simulation scheme was proposed to realize lattice gauge theories in 3+1 dimensions including dynamical fermions using only two-body interactions. Its main feature is the ability of obtaining the magnetic plaquette interactions without using fourth-order perturbation theory, thus resulting in stronger interactions and allowing the study of wider phase-space regions compared to analogue approaches. Second, following the aforementioned simulation scheme, an implementation of a lattice gauge theory with a non-abelian gauge groupthe dihedral group D 3 -was proposed, using ultracold atoms in optical lattices. Since the time evolution is performed in a Trotterized manner, intrinsic errors occur. These were studied in detail as a good bound on the trotter error gives more leeway to experimental errors.
The key ingredient of the digital simulation scheme is an auxiliary system which can be entangled with the physical system. This allows to create an isometry which mediates the complicated three and four-body terms of lattice gauge theory via the auxiliary system by using two-body interactions, as desired for implementations with various quantum simulation platforms. Moreover, it should be emphasized that all time evolutions in this algorithm are individually gauge invariant. The corresponding gauge group has to be either a compact Lie group or a finite group which is not a restriction for all relevant theories. In the case of compact Lie groups, the local Hilbert spaces of the gauge fields have an infinite dimension and therefore need to be truncated for a feasible implementation. However, since the isometry is defined in terms of the group element basis, the truncation has to be done there as well and can not be done in the typically used representation basis (see Sec. 2). Examples for such truncations are Z N for U (1) or -as proposed in this work -D N for O (2).
For the implementation of the lattice gauge theory with dihedral group D 3 -isomorphic to the symmetric group S 3 -we exploited the group structure of D 3 as a semidirect product. This allowed us to represent the gauge fields by a tensor product of a threelevel and a two-level system and thus simplified the implementation. The potential gain from this procedure would be even higher for more complicated gauge groups exhibiting a semidirect product structure.
No sophisticated experimental techniques (e.g. Feshbach resonances) are required. However, precise control over atomic collisions is needed in order to obtain the desired time evolution, in particular gates entangling the auxiliary system with the physical system, as they do not depend on the simulated time and are thus more prone to experimental errors. Future efforts on experimental techniques can therefore be targeted at the controllability of the relevant parameters, i.e. in particular fine tuning of the overlap integrals and the interaction time during scattering processes. The generation and experimental control of superlattices is important as well in order to create a staggering potential for the dynamical fermions. Also conducting experiments on simpler models -as currently set up for the Schwinger model -is a promising direction as it can serve as a proof of principle for the validity of quantum simulations of lattice gauge theories and might encourage more work in this direction. From the theoretical point of view, a logical next step is to think of possibilities to realize more complicated gauge groups. One step towards that goal is to find suitable ways to truncate compact gauge groups like for example SU (2) in a meaningful manner. The second term can be viewed as the electric energy of Z N as it acts trivially on the gauge field Hilbert space corresponding to reflections.

Appendix B. Trotter errors
For the bounds on the trotter error of the digital quantum simulation (presented in Sec. 3) a computation of commutators and nested commutators of the different parts of the Hamiltonian is required. Since the calculation of these commutators for a general lattice gauge theory is very lengthy, it is only sketched here.
Appendix B.0.1. First order For the first order formula the ordinary commutators need to be computed. Starting with the commutator between gauge-matter interactions on different links i and j, we obtain: (B.1) where we used the unitary operators U W from Sec. 3.2.2 to reduce the gauge-matter terms to pure fermionic tunneling terms, thus allowing to estimate this expression: where d U is the dimension of the representation of U under the gauge group and therefore the operator norm of the tunneling term. In the next step, the commutator between the matter-and gauge-matter interactions is calculated: The last commutator is the one between the magnetic and electric Hamiltonian. Since every link is contained in 2(d − 1) plaquettes, the commutator is straightforwardly estimated as: