Digital quantum simulation of lattice gauge theories in three spatial dimensions
Julian Bender, Erez Zohar, Alessandro Farace, J. Ignacio Cirac
DDigital quantum simulation of lattice gauge theoriesin three spatial dimensions
Julian Bender , Erez Zohar , Alessandro Farace and J. IgnacioCirac Max-Planck-Institut f¨ur Quantenoptik, Hans-Kopfermann-Straße 1, 85748 Garching,Germany.
Abstract.
In the present work, we propose a scheme for digital formulation oflattice gauge theories with dynamical fermions in 3+1 dimensions. All interactions areobtained 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 withoutthe use of perturbation theory, thus resulting in stronger interactions compared withanalogue approaches. The simulation scheme is applicable to lattice gauge theorieswith either compact or finite gauge groups. The required bounds on the digitizationerrors in lattice gauge theories, due to the sequential nature of the stroboscopic timeevolution, are provided. Furthermore, an implementation of a lattice gauge theorywith a non-abelian gauge group, the dihedral group D , is proposed employing theaforementioned simulation scheme using ultracold atoms in optical lattices.
1. Introduction
Gauge theories lie at the core of fundamental physics; the standard model of particlephysics - describing electromagnetic, weak and strong interactions - is based on theprinciple of gauge invariance [1]. It requires introducing additional degrees of freedom,the gauge fields, to the matter fields: force carriers, mediating interactions betweenmatter particles. If the coupling is small enough, perturbative expansions allowcalculations up to arbitrary accuracy, as in QED (Quantum Electrodynamics). In somequantum 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 (QuantumChromodynamics) at low energies. In such non-perturbative regimes only specialmethods can produce meaningful results.The most common approach so far has been lattice gauge theory [4, 5]. The idea isto discretize space (or spacetime) to construct a framework in which numerical toolscould be applied - with Monte Carlo methods being the most prominent ones [6]. Inspite 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 gaugetheories. A major one is the sign problem , which prevents investigations in fermionic a r X i v : . [ qu a n t - ph ] A p r igital quantum simulation of lattice gauge theories in three spatial dimensions igital quantum simulation of lattice gauge theories in three spatial dimensions e − itH = lim N →∞ ( (cid:81) j e − itH j /N ) N [51]. This method has already been proposed in 2+1dimensions to construct a digital scheme for lattice gauge theories with arbitrary gaugegroups [31]. A concrete quantum simulation with ultracold atoms has been proposedfor the groups Z and Z [31, 32].In this work we extend this proposal of an algorithm digitizing lattice gauge theorieswith arbitrary gauge groups to 3+1 dimensions. This is an important step towardsthe simulation of phenomena occurring in nature. To study the accuracy of the digitalscheme, a thorough analysis of the digitization (Trotter) error is conducted. Anotherimportant goal is the simulation of gauge theories with non-abelian gauge groups. Thesecond part of this work is therefore devoted to an ultracold atom based implementationof a lattice gauge theory with the simple non-abelian gauge group D , following thegeneral algorithm presented in the first part.This paper is organized as follows: First, a brief lattice gauge theory background will beprovided, with an emphasis on the Hamiltonian formulation used later on for quantumsimulation. In the second section the digital algorithm enabling quantum simulation oflattice gauge theories with dynamical fermions in three dimensions will be described.Afterwards, improved bounds on the digitization errors in lattice gauge theories will begiven. In the last section, possible implementations based on ultracold atoms will bediscussed, in particular the implementation of a lattice gauge theory with the dihedralgauge group D , by exploiting its semidirect product structure.
2. Hamiltonian formulation of lattice gauge theories
Lattice gauge theories can be formulated in a Hamiltonian framework exhibiting acontinuous time coordinate, as first proposed by Kogut and Susskind [52]. The latticeconsists 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 ∈ { , .., d } denotes the directionin which the link points). igital quantum simulation of lattice gauge theories in three spatial dimensions G to be either compact or finite, we label its irreduciblerepresentations by j and represent the matter fields by spinors ψ † jm , where m denotesthe components of j . Their behavior under group transformations, implemented by theunitary operator θ g , is (summing over repeated indices): θ g ψ † jm θ † g = ψ † jn D jnm ( g ) (1)where D jnm ( g ) is the irreducible unitary representation j of the group element g . We willwork with staggered fermions [53], distributing the Lorentz components of the spinorover neighboring lattice sites such that occupied even sites will correspond to particlesand vacant odd sites to anti-particles. The Dirac spinor is then regained in a continuumlimit. The gauge transformations ˇ θ g of staggered fermions are related to θ g byˇ θ g ( x ) = (cid:40) θ g for x ∈ eθ g det( D ( g − )) for x ∈ o (2)with e (o) denoting the even (odd) sublattice. We can define a state | D (cid:105) invariant underthe above transformation (analogous to the Dirac sea in the continuum) where all oddsites are fully occupied and all even sites are vacant.The other physical ingredients, the gauge degrees of freedom, are described by a tensorproduct of local Hilbert spaces on the links. The elements of each single link Hilbertspace can be expressed in the group element states {| g (cid:105)} g ∈ G . The group G can act onit in two ways, corresponding to left ( L ) and right ( R ) transformations:Θ Lg | h (cid:105) = | g − h (cid:105) , Θ Rg | h (cid:105) = | hg − (cid:105) (3)We define the group element operator U , a matrix of operators acting on the link Hilbertspace: U jmn = (cid:90) D jmn ( g ) | g (cid:105) (cid:104) g | d g (4)where for continuous groups d g is understood as the group (Haar) measure, whereas fordiscrete 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( U jmn ) † = (cid:90) dg | g (cid:105) (cid:104) g | ¯ D jmn ( g ) = (cid:90) dg | g (cid:105) (cid:104) g | D j † nm ( g ) = ( U j † ) nm (5)The group element operators obey the following rules under group transformations:Θ Lg U jmn Θ Lg † = D jmm (cid:48) ( g ) U jm (cid:48) n , Θ Rg U jmn Θ Rg † = U jmn (cid:48) D jn (cid:48) n ( g ) (6)(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 igital quantum simulation of lattice gauge theories in three spatial dimensions Lg ( x ,
2) Θ Lg ( x , Lg ( x ,
3) Θ Rg † ( x − , Rg † ( x − ,
1) Θ Rg † ( x − , θ † g ( x ) x + + + Figure 1.
The local gauge transformation Θ g ( x ), acting on the vertex x and adjacentlinks (shown here in three dimensions): ˇ θ † g ( x ) acts on the fermionic Fock space atvertex x , taking into account the staggered structure of the fermions. The three links( x , k ) emanating from vertex x are transformed by left transformations Θ Lg , whereasthe incoming links ( x − k , k ) are transformed by right transformations Θ Rg † . hand we can define a local gauge transformation which acts on all degrees of freedomintersecting at a vertex. It depends on a group element g which itself can depend onthe position (see Fig. 1 for illustration):Θ g ( x ) = (cid:89) k =1 ..d (cid:16) Θ Lg ( x , k )Θ Rg † ( x − k , k ) (cid:17) ˇ θ † g ( x ) (7)where k is the unit vector in k -direction. A state | ψ (cid:105) is therefore said to be gauge-invariant if Θ g ( x ) | ψ (cid:105) = | ψ (cid:105) , ∀ x , g (8)Introducing the dual basis to the group element states, the representationbasis {| jmn (cid:105)} , connected by the relation (cid:104) g | jmn (cid:105) = (cid:113) dim ( j ) | G | D mn ( g ) (with j labeling irreducible representations and m , n the components under left and righttransformations), we can define a gauge-invariant ”empty” state for the whole lattice,including matter and gauge fields: | (cid:105) ≡ | D (cid:105) (cid:79) links | (cid:105) (9)where | (cid:105) is a singlet state of the gauge fields in the representation basis, correspondingto the trivial representation. All other gauge invariant states can be obtained by actingwith gauge invariant operators on this trivial state. A conventional lattice gauge theoryHamiltonian consists of four such types of terms: igital quantum simulation of lattice gauge theories in three spatial dimensions magnetic Hamiltonian One can obtain gauge invariant operators by taking products of U -operators alongclosed paths. The shortest such possible path is a plaquette, characterized by twodirections k and l ( k < l and l ∈ { , .., d } ). Adding over all pairs of k and l forevery vertex x , one may construct: H B = λ B (cid:88) x ,k Introducing staggered fermions gives rise to the following staggered mass term: H M = M (cid:88) x ( − x ψ † n ( x ) ψ n ( x ) (13)where the alternating sign comes from the Dirac sea picture: particles on even sitesand anti-particles on odd sites.(iv) The gauge-matter Hamiltonian The last term is a fermionic hopping term minimally coupled to the gauge fields ina gauge invariant way: H GM = λ GM (cid:88) x ,k ψ † m ( x ) U mn ( x , k ) ψ n ( x + k ) + H.c. (14)The total Hamiltonian we want to simulate in the following chapters is the sum of allfour pieces. The state defined in (9) is the non-interacting vacuum: the ground state of H E + H M . igital quantum simulation of lattice gauge theories in three spatial dimensions 3. Digital algorithm for the quantum simulation of lattice gauge theories inthree dimensions Interactions in typical quantum simulation platforms are usually two-body, e.g. atomiccollisions in ultracold atomic setups or spin-spin interactions in trapped ion setups.Three-and four body processes are strongly suppressed on the relevant experimentaltimescales, making it much harder to map the Hamiltonian of the simulated model ontothe simulating system, if the former includes interactions of more than two bodies. Thisis particularly relevant for lattice gauge theories since magnetic interactions are four-body terms (see Sec. 2). For the purpose of quantum simulation of lattice gauge theoriesit is therefore desirable to design a scheme in which interactions involving three and moreconstituents can be rewritten as exact sequences of only two-body interactions. In thisway, the energy scale associated to plaquette interactions is not limited by perturbativearguments (as in previous proposals) and the simulation can give access to a biggerregion of the phase diagram.One approach to this problem is based on the idea of using an auxiliary degree offreedom that gets entangled with the physical degrees of freedom and mediates theirinteractions. In the following, we will briefly present an isometry which formalizes thisidea (it is sometimes referred to as stator [55, 56]). We anticipate that in this newframework the time evolution has to be realized stroboscopically due to the sequentialnature of the entangling operations with the auxiliary system. Therefore, a digitalalgorithm based on Trotter’s formula will be designed to simulate lattice gauge theoriesin three spatial dimensions, using only two-body interactions. This corresponds to thehybrid simulation scheme discussed in the introduction, where the time evolution isTrotterized but the individual parts of the Hamiltonian are still implemented by ananalogue simulation. In the last section, bounds on the Trotter error will be provided. 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 auxiliarydegrees of freedom (sometimes called control in the following). We denote the operatorsacting on the Hilbert space H by O ( H ). An isometry S can then be defined, mapping H A → H A ⊗ H B , which can be created by a unitary U AB ∈ O ( H A ⊗ H B ) acting on someinitial state | in B (cid:105) ∈ H B : S = U AB | in B (cid:105) ∈ O ( H A ) ⊗ H B (15)This can be viewed as an entangling operation between the physical and the auxiliarydegrees of freedom. If this entangling procedure is chosen in a certain way, operationson 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 (cid:48) ∈ H B in theauxiliary Hilbert space in such a way that the following relation holds: igital quantum simulation of lattice gauge theories in three spatial dimensions H (cid:48) S = SH (16)An analogue relation for the time evolution follows directly, since H (cid:48) n S = SH n : e − iH (cid:48) t S = Se − iHt (17)Therefore, by creating such an isometry and acting with H (cid:48) on the control, we obtainthe desired time evolution of the physical state | ψ A (cid:105) : e − iH (cid:48) t U AB | ψ A (cid:105) | in B (cid:105) = e − iH (cid:48) t S | ψ A (cid:105) = Se − iHt | ψ A (cid:105) (18)The evolved physical state is still entangled with the auxiliary state which means thatone can either perform another operation using the isometry S or disentangle bothstates. This would lead to a product state with the auxiliary state going back to itsinitial state: U † AB e − iH (cid:48) t U AB ( | ψ A (cid:105) ⊗ | in B (cid:105) ) = | in B (cid:105) ⊗ e − iHt | ψ A (cid:105) (19) In this section we discuss an algorithm to simulate the lattice gauge theory Hamiltonianin three spatial dimensions. We start from the lattice model described in Sec. 2. Tocreate plaquette and gauge-matter interactions by means of isometries, we introduce anauxiliary degree of freedom in the middle of every second cube (either all even or oddones) and assign to it a Hilbert space (cid:101) H isomorphic to the Hilbert spaces on the links(see Fig 2). Then, the lattice gauge theory Hamiltonian is split up into several partswhich are implemented independently and sequentially: H LGT = H E + H M + (cid:88) i =1 H B,i + (cid:88) j =1 H GM,j (20)where we explicitly distinguish gauge-matter interactions taking place along differentdirections and in odd or even cubes, as well as plaquette interactions corresponding tothe different plaquettes of a unit cube (therefore we get a sum of six terms in bothcases). The desired time evolution e − itH LGT is then approximated by a Trotterized timeevolution consisting of N steps: e − itH LGT ∼ ( (cid:81) j e − itH j /N ) N , where H j is any of the termsappearing in (20). While electric and mass terms can be treated easily using only thephysical degrees of freedom, the plaquette and gauge-matter terms are further decom-posed as a suitably chosen sequence of simpler interactions mediated by the auxiliarysystems. This sequence will then be executed in parallel for all cubes where auxiliarydegrees of freedom are located. However, since for the gauge-matter interactions theindividual parts of this sequence do not commute for adjacent links, we have to placethe auxiliary d.o.f. in every second cube to avoid undesired interactions. The exact igital quantum simulation of lattice gauge theories in three spatial dimensions Figure 2. The physical system consists of the gauge fields residing on the links (blue)and the matter fields on the vertices (red). The auxiliary degrees of freedom (green)are located in the center of every second cube (either even or odd). decompositions will be given in the next sections. Since we put auxiliary atoms in every second cube, wecan not realize all plaquette interactions at once and we split them up in the followingway: H B = (cid:88) x (cid:0) λ B Tr( U ( x , U ( x + , U † ( x + , U † ( x , H.c. (cid:1) + (cid:0) λ B Tr( U ( x , U ( x + , U † ( x + , U † ( x , H.c. (cid:1) + (cid:0) λ B Tr( U ( x , U ( x + , U † ( x + , U † ( x , H.c. (cid:1) ≡ (cid:88) x ( H B, ( x ) + H B, ( x ) + H B, ( x ))= (cid:88) x even ( H B, ( x ) + H B, ( x ) + H B, ( x )) + (cid:88) x odd ( H B, ( x ) + H B, ( x ) + H B, ( x )) ≡ H B, e + H B, e + H B, e + H B, o + H B, o + H B, o (21)It is important to mention that the six magnetic terms commute, therefore e − iτH B = (cid:81) j e − iτH B,je e − iτH B,jo and this splitting does not affect the error of the Trotterapproximation (20). To implement each term we will use the isometry S i = (cid:90) d g | g (cid:105) i (cid:104) g | i ⊗ | ˜ g (cid:105) (22) igital quantum simulation of lattice gauge theories in three spatial dimensions i , and thesecond one to the aforementioned auxiliary degree of freedom in the center of the cube.It fulfills the relation (cid:101) U S i = S i U link i (23)allowing to realize operations on the link i through the auxiliary degree of freedom. Theisometry S i can be created by the unitary U i = (cid:90) d g | g (cid:105) i (cid:104) g | i ⊗ Θ Lg † (24)acting on the initial state | (cid:101) in (cid:105) = | ˜ e (cid:105) . 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 S (cid:3) ( x ) = U (cid:3) ( x ) | (cid:101) in (cid:105) = U ( x ) U ( x ) U † ( x ) U † ( x ) | (cid:101) in (cid:105) (25)The crucial part is that it fulfills the relationTr( (cid:101) U ( x ) + (cid:101) U † ( x )) S (cid:3) ( x ) = S (cid:3) ( x ) Tr( U ( x ) U ( x ) U † ( x ) U † ( x ) + H.c. ) (26)i.e. acting locally with (cid:101) H B ( x ) = λ B Tr( (cid:101) U ( x ) + (cid:101) U † ( x )) (27)on the control of cube x enables us to realize the magnetic time evolution for thisplaquette. The required sequence acting on the plaquette state | ψ (cid:105) , the tensorproduct of the four link states, and the auxiliary state | (cid:101) in (cid:105) is U † (cid:3) ( x ) e − i (cid:101) H B ( x ) τ U (cid:3) ( x ) | ψ (cid:105) | (cid:101) in (cid:105) = | (cid:101) in (cid:105) e − iH B, ( x ) τ | ψ (cid:105) (28)The other two plaquette terms associated to cube x can be created in the samemanner but with different isometries. Using the abbreviations for the gauge fieldoperators defined according to Fig. 3, we need to replace the isometry S (cid:3) ( x ) by S (cid:3) ( x ) = U (cid:3) ( x ) | (cid:101) in (cid:105) (green, dashed plaquette), or S (cid:3) ( x ) = U (cid:3) ( x ) | (cid:101) in (cid:105) (blue,dotted plaquette). Applying the sequence from (28) gives then rise to the time evolutionof the physical state under H B, ( x ), or H B, ( 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 followingthree steps:(i) Create the isometry: Let all the controls interact with all the gauge fields on linksof type 4 and create the unitary (cid:81) x even U † ( x ). Repeat similar processes with links3, 2 and 1 to obtain the unitaries (cid:81) x even U † ( x ), (cid:81) x even U ( x ), (cid:81) x even U ( x ). In total, weget: (cid:81) x even U ( x ) U ( x ) U † ( x ) U † ( x ) = (cid:81) x even U (cid:3) ( x ). igital quantum simulation of lattice gauge theories in three spatial dimensions U U U U U U U U U xx + + + U ( x ) ≡ U ( x , U ( x ) ≡ U ( x + , U ( x ) ≡ U ( x + , U ( x ) ≡ U ( x , U ( x ) ≡ U ( x , U ( x ) ≡ U ( x + , U ( x ) ≡ U ( x + , U ( x ) ≡ U ( x + , U ( x ) ≡ U ( x + , Figure 3. There are three different plaquette terms associated to every vertex x : H B, ( x ) (red, solid plaquette), H B, ( x ) (green, dashed plaquette) and H B, ( x ) (blue,dotted plaquette). Each term involves four gauge field operators U , abbreviated asabove for a convenient description. (ii) Act on the controls with the Hamiltonian (cid:80) x even (cid:101) H B ( x ) for time τ , resulting in thetime evolution (cid:81) x even e − i (cid:101) H B ( x ) τ .(iii) In the last step, undo the isometry by creating the inverse of the first step, i.e. (cid:81) x even U † (cid:3) ( x ).The above procedure is applied to a state | ψ (cid:105) | (cid:101) in (cid:105) . Thanks to relation (28) we obtain: (cid:89) x even U † (cid:3) ( x ) e − i (cid:101) H B ( x ) τ U (cid:3) ( x ) | ψ (cid:105) | (cid:101) in (cid:105) = | (cid:101) in (cid:105) e − iH B, e τ | ψ (cid:105) ≡ | (cid:101) in (cid:105) W B, e | ψ (cid:105) (29)We repeat the procedure with the two isometries S (cid:3) and S (cid:3) . In this waywe create W B, e = e − iH B, e τ and W B, e = e − iH B, e τ . The same steps are then re-peated with the auxiliary degrees of freedom moved to the center of the odd cubesso that we can implement W B, o , W B, o , W B, o . Since all pieces of the magneticHamiltonian commute, this sequence gives us exactly the magnetic time evolution: W B, e W B, e W B, e W B, o W B, o W B, o = W B = e − iτH B . After expressing the four-body plaquette interac-tions as a sequence of two-body interactions, we want to obtain the gauge-matter inter-actions in a similar way. We need again to split up the relevant Hamiltonian terms intoparts suitable for implementation: igital quantum simulation of lattice gauge theories in three spatial dimensions H GM = (cid:88) x (cid:88) k =1 λ GM ψ † m ( x ) U mn ( x , k ) ψ n ( x + k ) + H.c. = (cid:88) x even ( H GM ( x , 1) + H GM ( x , 2) + H GM ( x , (cid:88) x odd ( H GM ( x , 1) + H GM ( x , 2) + H GM ( x , ≡ H GM, e + H GM, e + H GM, e + H GM, o + H GM, o + H GM, o (30)An important ingredient for rewriting these interactions as two-body terms is thefollowing unitary operation, entangling the fermion at vertex x and the gauge fieldon link ( x , k ): U W ( x , k ) = e iZ mn ( x ,k ) ψ † m ( x ) ψ n ( x ) (31)where Z mn = − i (log mat ( U )) mn , and the logarithm is taken only in matrix space (well-defined since the matrix elements commute). Its meaning becomes more apparent if weassume the gauge group G to be compact; then, we obtain U W ( x , k ) = e i ˆ φ a ( x ,k ) ψ † m ( x ) T amn ψ n ( x ) = e i ˆ φ a Q a (32)an interaction of the ”vector potential” operator ˆ φ a with the fermionic charge Q a . It cantherefore be interpreted as a fermionic transformation whose parameter is an operatoracting on the gauge field. The idea is now to use this transformation to map a purefermionic tunneling term into the desired gauge-matter interactions, as U W ( x , k ) ψ † n ( x ) U † W ( x , k ) = ψ † m ( x ) U mn ( x , k ) (33)Thus, defining the fermionic tunneling Hamiltonian as H t ( x , k ) = λ GM ( ψ † m ( x ) ψ m ( x + k ) + H.c. ) (34)allows writing the Hamiltonian H GM as: H GM ( x , k ) = U W ( x , k ) H t ( x , k ) U † W ( x , k ) (35)Since every fermion is connected to six links in three dimensions we have to split upthe process in six steps as described in the beginning. We start by realizing H GM, e , 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 thelink to obtain the unitary: (cid:81) x even U † W ( x , 1) .(ii) Allow tunneling on these links for time τ : (cid:81) x even e − iH t ( x , τ .(iii) Let the link degrees interact again with the fermions to generate: (cid:81) x even U W ( x , 1) . igital quantum simulation of lattice gauge theories in three spatial dimensions U ml ( x , U mn ( x , U mk ( x , ψ † m ( x ) ψ † l ( x + ) ψ † n ( x + ) ψ † k ( x + ) Figure 4. There are three gauge-matter terms associated to every vertex x ,corresponding to the three links emanating from this vertex: H GM ( x , 1) (red, solidlink), H GM ( x , 2) (green, dashed link) and H GM ( x , 3) (blue, dotted link). Eachinteraction consists of two fermions ψ † m ( x ) and ψ † n ( x + k ) located at the endpointsof the links and the gauge field operator U mn ( x , k ) on the link. This gives us in total (cid:89) x even U W ( x , e − iH t ( x , τ U † W ( x , 1) = e − i (cid:80) x even H GM ( x , τ ≡ W GM, e (36)By applying a similar sequence for the other links of the cube, we can create W GM, e , W GM, e , W GM, o , W GM, o , W GM, o .Using isometries, there is an alternative way of realizing the gauge-matterinteractions. It requires more steps but on the other hand does not require interactionsbetween the physical degrees of freedom as all of them are mediated by the auxiliarydegrees of freedom. The sequence goes as follows:(i) Let the controls - initially placed in all even cubes in the state | (cid:101) in (cid:105) = | ˜ e (cid:105) - interactwith the gauge links U according to (24) to create the isometry S : (cid:81) x even U ( x ) .(ii) Let the control interact with the fermion at vertex x to realize the interaction (cid:81) x even (cid:101) U † W ( x , 1) which is the same interaction as U † W ( x , 1) but between the controland the fermion ψ m ( x ). Due to the properties of the isometry S the interactionbetween the control and the fermion will translate into an interaction between thefermion and the link.(iii) Afterwards, allow for pure tunneling between the fermions which gives rise to (cid:81) x even e − iH t ( x , τ . igital quantum simulation of lattice gauge theories in three spatial dimensions (cid:101) U W ( x , 1) for all even cubes which is again realized by aninteraction between the control and the fermion ψ m ( x ): (cid:81) x even (cid:101) U W ( x , 1) .(v) Finally, we have to undo the isometry between the control and the gauge field: (cid:81) x even U † ( x ) .The resulting sequence - applied to some physical state | ψ (cid:105) and the auxiliary state | (cid:101) in (cid:105) - is: (cid:89) x even U † ( x ) (cid:101) U W ( x , e − iH t ( x , τ (cid:101) U † W ( x , U ( x ) | ψ (cid:105) | (cid:101) in (cid:105) = | (cid:101) in (cid:105) e − iH GM, e | ψ (cid:105) = | (cid:101) in (cid:105) W GM, e | ψ (cid:105) (37)We repeat a similar procedure for all other links in the cube which gives us W GM, e , W GM, e , W GM, o , W GM, o , W GM, o . The electric part W E = e − iH E τ and the matterpart W M = e − iH M τ are local terms of our Hamiltonian and thus one can implementthem by acting locally on the physical degrees of freedom.We can now write down the whole sequence for a time step τ (combining commutingmagnetic terms to W B ): W τ = W M W E W GM, o W GM, o W GM, o W GM, e W GM, e W GM, e W B (38)It’s important to notice that all time evolutions in the above sequence are individuallygauge-invariant. Therefore, errors coming from the digitization do not break gauge-invariance. Although the approximated Trotter evolution has the correct gauge symmetry, it isstill 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 digitizationscheme presented in the previous section. We focus on the standard Trotter formula (thefirst order formula) and the second order formula which gives a better approximationwithout 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 ofthe experimental parameters would have to be much more flexible and the number ofoperations required for a single time step would increase exponentially with the orderof the approximation [57]. The first order formula [51] is of the form U N ( t ) = ( (cid:89) j e − iH j tN ) N (39)For which, using the operator norm, the difference to the physical time evolution U ( t ) = e − itH can be bounded by [58–60]: igital quantum simulation of lattice gauge theories in three spatial dimensions (cid:107)U ( t ) − U N ( t ) (cid:107) ≤ t N (cid:88) j By inspection of (40) we see that for an upper bound on thedigitization error of the standard trotter formula, the commutators among all differentparts of the Hamiltonian in (43) have to be evaluated, as well as their norms. Since the igital quantum simulation of lattice gauge theories in three spatial dimensions (cid:107)U ( t ) − U N ( t ) (cid:107)≤ t N (cid:32) (cid:107) [ H E , H B ] (cid:107) + (cid:107) [ H GM , H E ] (cid:107) + (cid:107) H GM , H M (cid:107) + d − (cid:88) k =1 2 d (cid:88) j = k +1 (cid:107) [ H GM,j , H GM,k ] (cid:107) (cid:33) = t d U N links N (cid:18) λ B λ E d − 1) max j | f ( j ) | + λ GM λ E max j | f ( j ) | + M λ GM + λ GM d − (cid:19) (44)where d is the number of spatial dimensions, d U the dimension of the representationof the group element operator U and N links the number of links in the lattice. Onemight think that operator norms involving H E are unbounded but, since we either workwith finite groups (whose number of irreducible representations is finite) or appropriatetruncations of infinite gauge field Hilbert spaces, the expression max j | f ( j ) | is finite, sothat we always obtain sensible error bounds. To bound the error of the second order formula we needto calculate nested commutators according to (42). Details on their calculation can befound in the Appendix. We provide here the final result: (cid:107)U ( t ) − U ,N ( t ) (cid:107)≤ t N links d U N (cid:20) λ E λ B max j | f ( j ) | ( d − (cid:18) λ E max j | f ( j ) | + λ B d U ( d − (cid:19) + λ GM λ E max j | f ( j ) | (cid:18) λ GM d U (2(2 d − 1) + 1) + λ E max j | f ( j ) | (cid:19) + λ GM M (4 dλ GM + M ) + λ GM (2 d − (cid:18) 13 (4 d − 1) + 12 (cid:19)(cid:21) (45)If we assume a cubic lattice with L lattice sites per side we can express the numberof links as: N links = d ( L − L d − . The upper bound shows that N should scale as N ∼ L d/ t which is somewhat bad since it considers a very general setting. If werestrict ourselves to the observation of intensive quantities we expect this scaling tobe much better. However, there are observables in lattice gauge theories, e.g. Wilsonloops, which do not fulfill this requirement and thus need to be bounded by more generalestimates like the ones given above. 4. Implementation of digital lattice gauge theories with ultracold atoms With this general scheme for the digital construction of three-dimensional lattice gaugetheories at hand, we can turn to the implementation of some concrete examples withultracold atoms. Typical gauge groups of interest are compact (e.g. U (1)), for which igital quantum simulation of lattice gauge theories in three spatial dimensions representation basis [13, 14]. This procedure, however, spoilsunitarity 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 whichconverge to U (1) in the N → ∞ limit. The digital quantum simulation of Z N gaugetheories has been studied in [31, 32]. We summarize below their main features, and thenwe build on these to tackle the simulation of simple non-abelian gauge models withdihedral symmetry given by the group D N . Z N Lattice gauge theories with a finite abelian gauge group play an important role as theyapproximate compact quantum electrodynamics. Since the Hilbert space of the gaugefield is reduced to dimension N if the gauge group Z N is considered, ultracold atomscan be used to represent these gauge degrees of freedom. These N states are labeled by | m (cid:105) and we define unitary operators P and Q on them: P N = Q N = 1 P QP † = e i πN QQ | m (cid:105) = | m + 1 (cid:105) (cyclically) P | m (cid:105) = e i πN m | m (cid:105) (46)Since the group is abelian, its representations are one dimensional and we need toconsider a single fermionic species, ψ † , on the vertices. We can now define theHamiltonian of Z N lattice gauge theory with fermionic matter: H E = λ E (cid:88) x ,k (cid:0) − P ( x , k ) − P † ( x , k ) (cid:1) H B = λ B (cid:88) x ,k 1) and m ∈ { , }} (50)The above notation already suggests that D N can be decomposed into a semidirectproduct of the abelian groups Z N and Z corresponding to rotations and reflections: D N (cid:39) Z N (cid:111) Z . It is thus useful to write the states of the gauge field Hilbert spaceas states living in the tensor product of an N -dimensional Hilbert space and a two-dimensional one, | p, m (cid:105) = | p (cid:105) ⊗ | m (cid:105) ∈ H N ⊗ H . In the implementation, such aproduct Hilbert space can be realized by using two atoms with the appropriate hyperfinestructure. If we choose to work with the smallest faithful irreducible representation ofthe group, we need two different fermionic components for the matter, denoted by ψ and ψ , due to the non-abelian nature. Accordingly, the gauge field operators U on thelinks become matrices of operators U = e i πN ˆ pσ z σ ˆ mx (ˆ p acts on H N and ˆ m on H ; σ x and igital quantum simulation of lattice gauge theories in three spatial dimensions σ z act in matrix space). This allows us to write down the Hamiltonians H B = λ B (cid:88) x ,k 1) and ˜ c † m F ( x )(with m F = − / , / m F and do not change the internal level of the atoms. This can be achieved by liftingthe degeneracy of the hyperfine multiplets such that transitions changing m F will costenergy. A possible way to do this is by introducing a uniform magnetic field which addsthe following correction to the Hamiltonian (Zeeman shift): igital quantum simulation of lattice gauge theories in three spatial dimensions Figure 5. The simulating system consists of one atomic species on the verticesrepresenting the matter (red) and two for both the gauge fields (blue) on the linksand the controls (green) located at the center of every second cube (projected into twodimensions for better visualization). The simulated degrees of freedom are encoded inthe hyperfine structure of the atoms, i.e. either an F = 1 or an F = 1 / Dirac sea in the continuum). The empty green circles indicate the need to movethe auxiliary atoms between even and odd cubes. H Z = µ B g F m F B (56)where µ B is the Bohr magneton and g F the hyperfine Lande factor. The energy splittinghas to be different for different atomic species to avoid resonant exchanges, thereforewe need to choose species with different Lande factors. Another possible approachto realize the different energy splittings is to address each species individually, forexample exploiting the AC Stark effect. Second, at some point we need to modulate theinteraction strengths depending on the internal level of the atoms. This can be achievedfor 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 reachdifferent equilibrium positions within the same potential well. By properly choosing theLande factors of the atomic species and tuning the magnetic field gradient one can thentailor the overlap of the atomic Wannier wave functions (and hence their interactionstrength) in an m F -dependent way.Below we discuss several details of the implementation scheme. igital quantum simulation of lattice gauge theories in three spatial dimensions Before starting the simulationwe should define the initial configuration of our simulating system. It is reached if alloptical potentials are sufficiently deep and separated such that no tunneling occurs andall atomic wave functions do not overlap. All minima of the auxiliary lattice are loadedwith one atom in the group element state corresponding to the identity group element,i.e | (cid:101) in (cid:105) = | ˜ e (cid:105) = | ˜0 , ˜0 (cid:105) . This means we have to prepare the state | (cid:101) F = 1; ˜ m F = 0 (cid:105) for thethree-level system and | (cid:101) F = 1 / 2; ˜ m F = 1 / (cid:105) for the two level system. The preparationof the atoms representing gauge and matter fields depends on the initial physical statewe want to simulate. All atoms must occupy the motional ground state with energy E (different for different atomic species). As mentioned in the previous section, we alsointroduce a uniform magnetic field (or an AC Stark effect) to lift the degeneracy of theground state manifolds and induce energy splittings ∆ E (again different for differentspecies) between the different m F components.We can define the non-interacting Hamiltonian H which will be present throughout thewhole implementation: H = (cid:88) x ( E ,mat + ∆ E mat ) ψ † ( x ) ψ ( x ) + ( E ,mat − ∆ E mat ) ψ † ( x ) ψ ( x )+ (cid:88) x ,k (cid:88) m F ( E ,a + ∆ E a m F ) a † m F ( x , k ) a m F ( x , k )+ (cid:88) x ,k (cid:88) m F ( E ,c + ∆ E c m F ) c † m F ( x , k ) c m F ( x , k )+ (cid:88) x (cid:88) m F ( E , ˜ a + ∆ E ˜ a m F ) ˜ a † m F ( x )˜ a m F ( x )+ (cid:88) x (cid:88) m F ( E , ˜ c + ∆ E ˜ c m F ) ˜ c † m F ( x )˜ c m F ( x ) (57)All parts of the digital simulations are added on top of H . To recover the desiredHamiltonian H of our D lattice gauge theory, we move to an interaction picture , i.e.we will work in a rotating frame with respect to H and make use of the rotating waveapproximation. The mass Hamiltonian in three dimensions takes theform H M = M (cid:88) x ( − x + x + x ψ † ( x ) ψ ( x ) (58)with ψ † ( x ) ψ ( x ) = ψ † ( x ) ψ ( x ) + ψ † ( x ) ψ ( x ). Thus, the corresponding time evolution W M = e − iH M τ for a time step τ can be implemented by smoothly modulating thesuperlattice trapping the fermions so that the energy of the even minima is increasedby an amount M even . This results in the Hamiltonian igital quantum simulation of lattice gauge theories in three spatial dimensions H (cid:48) M = M even (cid:88) x (1 + ( − x + x + x ) ψ † ( x ) ψ ( x ) (59)If we act with this Hamiltonian for time M even M τ , we obtain the desired unitary evolution W M , up to an irrelevant global phase. The creation of plaquette interactions and gauge-matterinteractions involves constructing the isometry S i (see Sec. 3), entangling auxiliaryatoms with the atoms on link i . If we want to create it from the auxiliary statecorresponding to the neutral element | (cid:101) in (cid:105) = | ˜ e (cid:105) , we have to apply U i = (cid:82) dg | g (cid:105) i (cid:104) g | i ⊗ Θ L † g .Specifying this equation to the gauge group D , we obtain the following interactionbetween the d.o.f. on link i and the ones of the control: U i | (cid:101) in (cid:105) = (cid:88) p (cid:88) m | p, m (cid:105) i (cid:104) p, m | i ⊗ Θ L † p,m | ˜0 , ˜0 (cid:105) = (cid:88) p (cid:88) m | p, m (cid:105) i (cid:104) p, m | i ⊗ (cid:101) Q p (cid:101) Q m | ˜0 , ˜0 (cid:105) = (cid:101) Q ˆ p ,i (cid:101) Q ˆ m ,i | ˜0 , ˜0 (cid:105) (60)where ˆ p = (cid:80) p p | p (cid:105) i (cid:104) p | i and ˆ m = (cid:80) m m | m (cid:105) i (cid:104) m | i . As defined previously, (cid:101) Q and (cid:101) Q arethe raising operators of the auxiliary atoms, i.e (cid:101) Q ˆ m ,i and (cid:101) Q ˆ p ,i raise the ˜ m F -values of theauxiliary atoms according to the m F -values of the atoms on link i . By choosing | ˜0 , ˜0 (cid:105) as the initial state of the auxiliary atoms, the creation of the isometry reduces to aninteraction between the three-level atom on the link and the auxiliary three-level atomin parallel with an interaction between the two-level atom on the link and the auxiliarytwo-level atom. These are the same interactions required for creating the isometry ofa Z lattice gauge theory [31], respectively a Z lattice gauge theory [32]. We cantherefore directly adopt the procedure from [31,32]. The idea is to write (cid:101) Q ˆ p ,i and (cid:101) Q ˆ m ,i asan interaction between the z-components of the hyperfine angular momentum operators (cid:101) F z, and F z, , respectively (cid:101) F z, and F z, : (cid:101) Q ˆ p ,i = (cid:101) V † U (cid:48) ,i (cid:101) V with U (cid:48) ,i = e i π (cid:101) F z, F z, (61)where (cid:101) V † is a local change of basis from the (cid:101) P -basis {| ˜ p (cid:105)} to its conjugate (cid:101) Q -basisand: (cid:101) Q ˆ m ,i = (cid:101) V † U (cid:48) ,i (cid:101) V with U (cid:48) ,i = e − i π ( i π ) (1 − F z, )(1 − (cid:101) F z, ) (62)where (cid:101) V † is mapping from the (cid:101) P -basis {| ˜ m (cid:105)} into the conjugate (cid:101) Q -basis. The basistransformations (cid:101) V and (cid:101) V are local operations on the auxiliary atoms that can beimplemented with optical/RF fields. The interactions between the z-components ofthe hyperfine angular momentum operator can be realized by introducing an energysplitting between the different m F -levels such that the two-body scattering term willcontain only m F preserving terms. The sequence to obtain U i is therefore: U i = (cid:101) V † (cid:101) V † U (cid:48) ,i U (cid:48) ,i (cid:101) V (cid:101) V (63) igital quantum simulation of lattice gauge theories in three spatial dimensions m F = 1 and ˜ m F = − (cid:101) F z, into − (cid:101) F z, . In the same way, the conjugate of the two-level system is created. Knowing how to construct the isometry, theimplementation of the plaquette interactions is straightforward. Since we have to splitthem in six different parts (see Sec. 3), we start with H B, e , the type 1 plaquettes ofthe even cubes, where the auxiliary atoms are placed in the standard configuration. Wefollow 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 tailoringthe interactions as described above (neglecting the basis transformations (cid:101) V for themoment). This can be done in parallel for the whole lattice: U (cid:48) ie = (cid:89) x even U (cid:48) ,i ( x ) U (cid:48) ,i ( x ) (64)The desired plaquette isometry is obtained by applying this procedure to all fourlinks and including overall basis transformations (cid:101) V ,all and (cid:101) V ,all : (cid:89) x even U ( x ) U ( x ) U † ( x ) U † ( x ) = (cid:101) V † ,all (cid:101) V † ,all U (cid:48) e U (cid:48) e U (cid:48)† e U (cid:48)† e (cid:101) V ,all (cid:101) V ,all (65)This operation, acting on the tensor product of | ˜0 , ˜0 (cid:105) and any state of the links,gives rise to the proper entangled state which maps plaquette interactions to localoperations on the control.(ii) The next step is a local operation on the auxiliary Hilbert space. We need toimplement e − i (cid:101) H B τ with (cid:101) H B being the control Hamiltonian (cid:101) H B = λ B Tr ( (cid:101) U + (cid:101) U † ).This requires an interaction between the two-level and the three-level system: (cid:101) H B = λ B Tr (cid:32)(cid:88) p (cid:88) m | ˜ p, ˜ m (cid:105) (cid:104) ˜ p, ˜ m | e i π σ z p σ mx + H.c. (cid:33) = 2 λ B ( (cid:101) P + (cid:101) P † )(1 − ˆ m ) (66)where ˆ m = (cid:80) m m | ˜ m (cid:105) (cid:104) ˜ m | . We can rewrite (cid:101) P + (cid:101) P † = − I + 3 | ˜0 (cid:105) (cid:104) ˜0 | ≡ − I + 3 N with (cid:101) N ≡ ˜ a † ˜ a . Defining a number operator for the | (cid:101) F = 1 / m F = 1 / (cid:105) state ofthe two-level system as (cid:101) N / ≡ ˜ c † / ˜ c / we can write down the interaction e − i (cid:101) H B τ : e − i (cid:101) H B τ = e − i λ B ( − I +3 (cid:101) N ) (cid:101) N / τ = e i λ B (cid:101) N / τ e − i λ B (cid:101) N (cid:101) N / τ (67)The first exponential is a local term of the two-level system which can beimplemented by means of optical/RF fields. The second term requires scatteringbetween the two auxiliary atoms. The corresponding Hamiltonian density in secondquantized form is [63]: igital quantum simulation of lattice gauge theories in three spatial dimensions H scat ( x ) = 2 πµ (cid:88) α,β,γ,δ n − (cid:88) k =0 g k (( F · F ) k ) α,β,γ,δ Φ † α ( x )Φ † β ( x )Φ γ ( x )Φ δ ( x ) (68)where Φ † α denotes the creation operator of the atomic Wannier wave functioncorresponding to the internal state α and µ the reduced mass of the two atomicspecies. The projection operators onto the different scattering channels areexpressed by polynomials of F · F , the coefficients g k are therefore functions ofthe scattering lengths. To obtain the time evolution due to this interaction we haveto integrate the Hamiltonian density over space and time. Since eq. (68) involvesonly specific levels, we need to turn on the magnetic field gradient and split thedifferent m F components such that only the ˜ m F = 0-component and the ˜ m F = 1 / U scat, = I + ( e − ig α − | ˜0 , ˜ (cid:105) (cid:104) ˜0 , ˜ | = e − ig α (cid:101) N (cid:101) N / (69)with g = (3 a / + 4 a / ) ( a / , a / are the scattering lengths for the scatteringchannels with F tot = 1 / F tot = 3 / 2) and α the time-integral of the overlapof the two wave-functions during the collision. By carefully tuning the interactiontime we can set α = λ B τg and finally obtain: U scat, = e − i λ B (cid:101) N (cid:101) N / τ (70)which is up to local operations the desired unitary V B . This interaction will beimplemented in parallel for all cubes where auxiliary atoms are placed, i.e. in thiscase for the even cubes. Hence, the overall interaction of this step is e − i (cid:80) x even (cid:101) H B ( x ) τ .When the magnetic field gradient is on, different levels of the hyperfine multipletwill 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 andaccumulate phases in the opposite direction until the net effect is zero (this trickhas 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 takingthe hermitian conjugate of the first step, i.e. the sequence: (cid:101) V † ,all (cid:101) V † ,all U (cid:48) e U (cid:48) e U (cid:48)† e U (cid:48)† e (cid:101) V ,all (cid:101) V ,all (71)According to (29) these three steps give us W B, e . If we repeat now the same procedurebut with the links corresponding to the second and third plaquette term, we obtain W B, e and W B, e . To realize the odd cubes time evolution, we move the auxiliary atomsto the centers of the odd cubes and repeat all of the above. This results in the time igital quantum simulation of lattice gauge theories in three spatial dimensions W B, o , W B, o and W B, o . Afterwards, the auxiliary atoms are brought back tothe centers of the even cubes. For the Gauge-matter interactions on a link ( x , k )we have to implement the Hamiltonian H GM ( x , k ) = λ GM ψ † a ( x ) U ab ( x , k ) ψ b ( x + k ) + H.c. = λ GM ψ † a ( x ) ( e i π σ z ˆ p ) ab ( σ ˆ mx ) bc ψ c ( x + k ) + H.c. = λ GM ψ † a ( x ) ( U p ) ab ( x , k )( U m ) bc ( x , k ) ψ c ( x + k ) + H.c. (72)with U p ≡ e i π σ z ˆ p and U m ≡ σ ˆ mx . We can use the product structure of U to implementthe 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 : U W,p ( x , k ) = e log( U p ) ab ( x ,k ) ψ † a ( x ) ψ b ( x ) = e i π ˆ p ( ψ † ( x ) ψ ( x ) − ψ † ( x ) ψ ( x )) (73)and another one corresponding to U m : U W,m ( x , k ) = e log( U m ) ab ( x ,k ) ψ † a ( x ) ψ b ( x ) = e i π ˆ m ( ψ † ( x ) ψ ( x )+ ψ † ( x ) ψ ( x ) − ψ † ( x ) ψ ( x ) − ψ † ( x ) ψ ( x )) (74)With these definitions at hand we can get the following relation by applying twice theBaker-Campbell-Hausdorff formula: U W,p ( x , k ) U W,m ( x , k ) ψ † n ( x ) U † W,m ( x , k ) U † W,p ( x , k ) = ψ † a ( x )( U p ) ab ( x , k )( U m ) bn ( x , k ) (75)The gauge-matter Hamiltonian can then be written as H GM ( x , k ) = U W,p ( x , k ) U W,m ( x , k ) H t ( x , k ) U † W,m ( x , k ) U † W,p ( x , k ) (76)with the tunneling Hamiltonian H t ( x , k ) = λ GM ( ψ † a ( x ) ψ a ( x + k ) + H.c. ) The crucialthing to note here is that all the terms involve only two-body interactions which allowsan implementation with the proposed ultracold atomic setup. We can not implementall gauge-matter interactions at once as the fermions on the vertices are only allowed tointeract with one link at a time. Focusing on the links in the -direction for the evencubes, we describe how to realize the time evolution e − i (cid:80) x even H GM ( x , τ . Since we want tokeep the lattice of the matter and link degrees of freedom fixed, these interactions willbe mediated by the control atoms according to the algorithm presented in 3.2.2.(i) We first build the isometry S between auxiliary atoms located at the center of evencubes x and the corresponding atoms on link ( x , (cid:81) x even U ( x ). This interactioncan be implemented exactly in the same way as already done for the plaquetteterm (see (64)). Due to the relation in (23) the gauge-matter interactions will thentranslate into an interaction of exactly the same form but between the auxiliaryatoms and the fermions. igital quantum simulation of lattice gauge theories in three spatial dimensions U † W,p and U † W,m have to be implemented by two-bodyscattering processes but between the fermions and the auxiliary atoms due to theisometry, therefore denoted as (cid:101) U W,p and (cid:101) U W,m . Starting with (cid:101) U † W,p , we first writeit in terms of the angular momentum operator respectively the second quantizedoperators ψ and ψ for the fermions: (cid:101) U † W,p ( x , k ) = e − i π ˆ p ( ψ † ( x ) ψ ( x ) − ψ † ( x ) ψ ( x )) = e − i π (cid:101) F z, ( ψ † ( x ) ψ ( x ) − ψ † ( x ) ψ ( x )) (77)Now we have to tailor the atomic collision between the (cid:101) F = 1 and the F = 1 / m F -values. The interactionHamiltonian contains two possible scattering channels and gives rise to the followingtime evolution: U scat, = e − iβ ( g ( ψ † ψ + ψ † ψ )+ g (cid:101) F z, ( ψ † ψ − ψ † ψ )) (78)with g = (3 a / + 4 a / ), g = ( a / − a / ) ( a / , a / are the scattering lengthsfor the scattering channels with F tot = 1 / F tot = 3 / 2) and β the time-integralof the wave-function overlap. If we tune overlap and interaction time such that β = π g we obtain U scat, = e − i πg g ( ψ † ψ + ψ † ψ ) e − i π (cid:101) F z, ( ψ † ψ − ψ † ψ ) (79)The second exponential is the desired interaction (cid:101) U † W,p whereas the first term is afermion-dependent phase, denoted from now on as V W (cid:48) ( θ ) = e − iθψ † ψ (80)where θ = πg g and ψ † ψ = ψ † ψ + ψ † ψ . A discussion of these phases will be donelater on. Before, the implementation of (cid:101) U † W,m is explained. It has the form: (cid:101) U † W,m = e − i π ˆ m ( ψ † ψ + ψ † ψ − ψ † ψ − ψ † ψ ) = e − i π (cid:101) N − / ( ψ † ψ + ψ † ψ − ψ † ψ − ψ † ψ ) = V H,fer e − iπ (cid:101) N − / ψ † ψ V H,fer (81)with (cid:101) N − / ≡ ˜ c †− / ˜ c − / and V H,fer = √ ( σ x,fer + σ z,fer ) a Hadamard transformon the fermions which can be implemented by means of optical/RF fields. Theremaining two-body interaction is realized as scattering between the F = 1 / U scat, = e − iγ ( g (cid:80) m ˜ c † m ˜ c m ( ψ † ψ + ψ † ψ )+ g (cid:101) F z, ( ψ † ψ − ψ † ψ ) ) (82)(for the explicit form of g k see [32]). We switch on a magnetic field gradientdesigned in a way that only the m F = − / γ = πg + g which gives rise to: U scat, = e − iγ (cid:16) g ˜ c †− / ˜ c − / ψ † ψ + g ˜ c †− / ˜ c − / ψ † ψ (cid:17) = e − iπ ˜ c †− / ˜ c − / ψ † ψ (83) igital quantum simulation of lattice gauge theories in three spatial dimensions (cid:101) U † W,p and (cid:101) U † W,m is done in parallel for all even cubes weget the sequence (cid:89) x even (cid:101) U † W,m ( x , (cid:101) U † W,p ( x , V W (cid:48) ( θ ) (84)(iii) In the next step we implement the tunneling in the 1-direction for even cubeswhich can be achieved by modulating the superlattice and decreasing the potentialbarriers on the desired links. We get (cid:89) x even e − iH t ( x , τ (85)(iv) After the tunneling we need to realize the conjugate of (cid:101) U † W,p and (cid:101) U † W,m , i.e. (cid:101) U W,p and (cid:101) U W,m . One way of creating (cid:101) U W,p is by doing a spin flipping operation (cid:101) V F, forthe three-level system of the control which results in: (cid:101) V F, (cid:101) U † W,p (cid:101) V † F, = (cid:101) U W,p (cid:101) V F, V W (cid:48) ( θ ) (cid:101) V † F, = V W (cid:48) ( θ ) (86)For the creation of (cid:101) U W,m we simply observe that (cid:101) U † W,m is hermitian. The sequencefor step 4 is (cid:89) x even V W (cid:48) ( θ ) (cid:101) U W,p ( x , (cid:101) U W,m ( x , 1) (87)(v) In the last step we need to undo the isometry, which is done by the conjugate ofthe first step, (cid:81) x even U † ( x ) (see Sec. 4.2.4).We summarize by writing down the whole sequence acting on the initial auxiliary state | (cid:101) in (cid:105) = | ˜0 , ˜0 (cid:105) : (cid:89) x even V W (cid:48) ( θ ) U † ( x ) (cid:101) U W,p ( x , (cid:101) U W,m ( x , e − iH t ( x , τ (cid:101) U † W,m ( x , (cid:101) U † W,p ( x , U ( x ) V W (cid:48) ( θ ) | (cid:101) in (cid:105) = | (cid:101) in (cid:105) V W (cid:48) ( θ ) (cid:89) x even e − iH GM ( x , τ V W (cid:48) ( θ ) = | (cid:101) in (cid:105) V W (cid:48) ( θ ) W GM, e V W (cid:48) ( θ ) (88)We finally get the desired gauge-matter interactions up to the fermionic phases V W (cid:48) ( θ ).However, if we consider the whole lattice (on which the number of fermions is globallyconserved) it can be shown that the phases correspond to a static vector potential ofzero magnetic field and are therefore unphysical, as carried out in the procedure givenin [31]. If we repeat the whole sequence (88) for the other links we obtain the gauge-matter interactions W GM, e , W GM, e , W GM, o , W GM, o and W GM, o . The electric Hamiltonian for the gauge group D acts onthe gauge fields residing on the links. If we choose its second part - which corresponds igital quantum simulation of lattice gauge theories in three spatial dimensions Z we obtain, usingthe notation of previous sections: H E = λ E (cid:88) x ,k h E ( x , k )with h E ( x , k ) = 12 f r (cid:88) m,m (cid:48) | , m (cid:105) (cid:104) , m (cid:48) | + f l (1 − P − P † ) ⊗ I (89)If we also express the interactions of the first part in terms of operators acting on thelink atoms, we end up with: h E ( x , k ) = 12 f r a † a ⊗ (1 + σ x ) + f l (cid:88) m F = − (1 + | m F | ) a † m F a m F ⊗ I (90)The first Hilbert space represents the three-level system, the second one the two-levelsystem. The coefficient f l is the overall coefficient for the electric part corresponding topure rotations, equivalently to Z . We have to implement the time evolution: W E = e − iH E τ = (cid:89) x ,k e − ih E ( x ,k ) τ (91)with e − ih E τ = e − i λEfr a † a τ e − i λEfr a † a σ x τ e − iλ E f l (cid:80) mF (1+ | m F | ) a † mF a mF τ (92)The first and the third exponential are local terms of the atoms and can be addressed byexternal fields. The second term is implemented by two-body scattering similar to theone 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 middleof the link. Following the steps for the plaquette interactions, we obtain: U scat, = e − iδg N N / (93)Tuning overlap and interaction time such that δ = λ E f r τg and combining it with the localoperation V = e i λEfrτ N / , gives us: V U scat, = e i λEfrτ N e − iλ E f r N N / τ = e − i λEfr N σ z τ (94)If we then perform a Hadamard transform V H, on the two-level system, we get thedesired interaction: V H, V U scat, V H, = e − i λEfr N σ x = e − i λEfr a † a σ x τ (95)which gives us the electric Hamiltonian up to local operations.We have implemented all interactions using local operations on the atoms and tailoringthe appropriate two-body scattering terms. If we use the sequence to evolve thesystem for a time τ = T /N and we repeat the same sequence N times, we get aTrotter approximation of the desired time-evolution e − iH LGT T . The accuracy of thisapproximation is discussed below. igital quantum simulation of lattice gauge theories in three spatial dimensions The errors affecting the precision of the simulation are twofold. On theone hand, we have Trotter errors coming from the digitization which can be estimatedby specifying the general error bounds given in Sec. 3 to the case of three dimensionand gauge group D . We obtain for the fist order formula (see (44)): (cid:107)U ( t ) − U N ( t ) (cid:107) ≤ t ( L − L N (cid:18) λ B λ E + 2 λ GM λ E + M λ GM + λ GM (cid:19) (96)and the second order formula (see (45)): (cid:107)U ( t ) − U ,N ( t ) (cid:107) ≤ t ( L − L N (128 λ E λ B ( λ E + λ B ) + 2 λ GM λ E (22 λ GM + λ E )+ λ GM M (6 λ GM + 12 M ) + 12512 λ GM (cid:19) (97)We stress again that the digitization error doesn’t break gauge invariance, because allsteps of the sequence individually respect the right symmetry. Therefore, the Trotterexpansion can only give rise to quantitative deviations, but not to qualitative changes.On the other hand, there will be experimental errors in the implementation. Unlikeerrors caused by the Trotterization, they may break the gauge symmetry and accumulatestep by step. We briefly want to look at the scaling of these errors. We considera small perturbation h j to one of the Hamiltonians H j which is realized during theimplementation 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 normby (cid:107) h j (cid:107) τ . 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 secondorder formula (41), i.e. the gate is repeated 2 N times). We have to distinguish fourcases: 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 andobtain a scaling of √ N with the number of Trotter steps compared to a linear scalingin the case of a systematic error. In the same vein, a gate depending on the simulatedtime t is advantageous since the time step τ = t N in each trotter sequence scales as 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 aresummarized in Table 1. We see that operations that do not depend on the simulatedtime t are the ones most prone to errors. During their implementation a lot of careshould be taken, in particular to avoid systematic errors. When estimating the errorof the whole implementation sequence, one should keep in mind that errors of differentgates are generally independent and thus do not add up linearly. However, the totalexperimental error will still increase with N , so that the number of Trotter steps has tobe chosen in a way to balance digitization and implementation errors.Typical sources of errors in ultracold atom experiments are as follows: The first oneis decoherence, e.g. caused by spontaneous scattering of lattice photons with the atoms, igital quantum simulation of lattice gauge theories in three spatial dimensions t (cid:107) h j (cid:107) t √ N Systematic error/dependence on t (cid:107) h j (cid:107) t Statistical error/no dependence on t (cid:107) h j (cid:107)√ N t exp ,j Systematic error/no dependence on t (cid:107) h j (cid:107) N t exp ,j Table 1. The different types of experimental errors corresponding to some gate j are distinguished by the nature of the perturbation h j (statistical or systematic) andwhether the gate depends on the simulated time t or is a fixed operation lasting forsome time t exp ,j . The error bound for each type scales differently with the number ofTrotter steps N . 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 orderof minutes have already been achieved [15, 64, 65], thus requiring the total simulationtime t sim to fulfill t sim << t coh . Secondly, one needs to ensure that the atoms remainin the lowest Bloch band throughout the whole implementation. Hence, it is of crucialimportance to shape the lattice and move the atoms in an adiabatic way. This isparticularly important in our simulation scheme, where the auxiliary atoms have to bemoved around or when the matter lattice has to be deformed to allow tunneling. Thismeans that the corresponding timescale t mov should be bigger than the inverse of thefrequency ω associated to the energy difference between lowest and first excited Blochband ( t mov >> /ω ), while at the same time the obvious constraint t mov << t coh has tobe fulfilled. However, such techniques have also become well-controlled [15, 66].Errors more specific to this proposal are connected with the tailoring of the two-body scattering. This requires a high degree of control over the overlap of the atomicwave functions and accurate timing of interaction during these collisions. This is alsodependent on the ability to design and manipulate the magnetic field gradient in aprecise manner. 5. Summary In this work, two main results were discussed. First, a digital simulation schemewas proposed to realize lattice gauge theories in 3+1 dimensions including dynamicalfermions using only two-body interactions. Its main feature is the ability of obtainingthe magnetic plaquette interactions without using fourth-order perturbation theory, thusresulting in stronger interactions and allowing the study of wider phase-space regionscompared to analogue approaches. Second, following the aforementioned simulationscheme, an implementation of a lattice gauge theory with a non-abelian gauge group -the dihedral group D - was proposed, using ultracold atoms in optical lattices. Since thetime evolution is performed in a Trotterized manner, intrinsic errors occur. These werestudied in detail as a good bound on the trotter error gives more leeway to experimentalerrors. igital quantum simulation of lattice gauge theories in three spatial dimensions 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 - isomorphicto the symmetric group S - we exploited the group structure of D as a semidirectproduct. This allowed us to represent the gauge fields by a tensor product of a three-level and a two-level system and thus simplified the implementation. The potential gainfrom this procedure would be even higher for more complicated gauge groups exhibitinga 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 desiredtime evolution, in particular gates entangling the auxiliary system with the physicalsystem, as they do not depend on the simulated time and are thus more prone toexperimental errors.Future efforts on experimental techniques can therefore be targeted at the controllabilityof the relevant parameters, i.e. in particular fine tuning of the overlap integrals and theinteraction time during scattering processes. The generation and experimental controlof superlattices is important as well in order to create a staggering potential for thedynamical fermions. Also conducting experiments on simpler models - as currently setup for the Schwinger model - is a promising direction as it can serve as a proof ofprinciple for the validity of quantum simulations of lattice gauge theories and mightencourage more work in this direction.From the theoretical point of view, a logical next step is to think of possibilities torealize more complicated gauge groups. One step towards that goal is to find suitableways to truncate compact gauge groups like for example SU (2) in a meaningful manner. Acknowledgments J.I.C. is supported by the ERC QENOCOBA under the EU Horizon 2020 program(grant agreement 742102). igital quantum simulation of lattice gauge theories in three spatial dimensions Appendix A. Details on D N lattice gauge theory In the following we will present some details on the lattice gauge theory of D N . Dueto its non-abelian gauge group the representations of the group become non-trivial andthus a lot of terms more complicated. Therefore, we will start by discussing the mostimportant group properties and the irreducible representations of D N . D N is the symmetry group of rotations by πN in a two-dimensional plane and reflectionsalong a certain axis (any axis passing through the center of rotations is possible). It canbe characterized by the set D N = { g = ( p, m ) ≡ R (2 π/N ) p S m | p ∈ [0 , N − 1) and m ∈ { , }} (A.1)The structure of the group is defined by the composition rules:( p, m ) · ( r, n ) = ( p + ( − m r, m + n ) (A.2)where the addition of p and r is understood as modulo N , respectively modulo 2 for m and n . The neutral element is e = (0 , 0) and the inverse element of ( p, m ) is( p, m ) − = ( p ( − m +1 , m ). The representation theory of D N ( N odd and N ≥ 3) ischaracterized by the three irreducible representations shown in the table below:Trivial (dimension 1) D t ( p, m ) = 1Sign (dimension 1) D s ( p, m ) = ( − m k-th (dimension 2) D k ( p, m ) = e i πpN kσ z σ mx We exclude the cases where N is even, since they have additional sign representationsand are not relevant for the discussion of D . With the above table, the electricHamiltonian can easily be given in the representation basis. However, since this form ofthe Hamiltonian is not very feasible for the proposed quantum simulation we will showhow to transform it to the group element states: h E ( x , k ) = (cid:88) g,g (cid:48) (cid:88) j,m,n f ( j ) | g (cid:105) (cid:104) g | jmn (cid:105) (cid:104) jmn | g (cid:48) (cid:105) (cid:104) g (cid:48) | = (cid:88) g,g (cid:48) (cid:88) j dim ( j ) | G | f ( j ) | g (cid:105) Tr( D j ( g ) D j † ( g (cid:48) )) (cid:104) g (cid:48) | (A.3)To specify this expression for D N we need to calculate the trace from above for all irre-ducible representations:Trivial representation: Tr( D t ( p, m ) D t † ( p (cid:48) , m (cid:48) )) = 1Sign representation: Tr( D s ( p, m ) D s † ( p (cid:48) , m (cid:48) )) = ( − m + m (cid:48) k-th representation: Tr( D k ( p, m ) D k † ( p (cid:48) , m (cid:48) )) = δ mm (cid:48) ( e i πN k ( p − p (cid:48) ) + e − i πN k ( p − p (cid:48) ) )Inserting this into (A.3) we obtain: igital quantum simulation of lattice gauge theories in three spatial dimensions h E ( x , k ) = 12 N (cid:88) p,p (cid:48) (cid:88) m,m (cid:48) | p, m (cid:105) (cid:104) p (cid:48) , m (cid:48) | (cid:32) f t + f s ( − m + m (cid:48) + 2 N − (cid:88) k =1 f k δ mm (cid:48) ( e i πN k ( p − p (cid:48) ) + e − i πN k ( p − p (cid:48) ) ) (cid:33) (A.4)The expression simplifies if we go to the conjugate basis of {| p (cid:105)} which can be viewedas the angular momentum basis {| l (cid:105)} characterized by the relation (cid:104) l, m | p, n (cid:105) = 1 √ N δ mn e − i πN lp . (A.5)We obtain h E ( x , k ) = 12 N (cid:88) m,m (cid:48) (cid:88) l,l (cid:48) | l, m (cid:105) (cid:88) p,p (cid:48) (cid:16) f t e − i πN lp e i πN l (cid:48) p (cid:48) + f s e − i πN lp e i πN l (cid:48) p (cid:48) ( − m + m (cid:48) +2 N − (cid:88) k =1 f k δ mm (cid:48) ( e i πN ( k − l ) p e i πN ( l (cid:48) − k ) p (cid:48) + e − i πN ( k + l ) p e i πN ( k + l (cid:48) ) p (cid:48) ) (cid:33) (cid:104) l (cid:48) , m (cid:48) | = 12 (cid:88) m,m (cid:48) (cid:32) f t | , m (cid:105) (cid:104) , m (cid:48) | + f s ( − m + m (cid:48) | , m (cid:105) (cid:104) , m (cid:48) | + 2 (cid:88) l (cid:54) =0 f l | l, m (cid:105) (cid:104) l, m | (cid:33) (A.6)where the coefficients f l have to satisfy the constraint f l = f − l ∀ l . If we redefine thecoefficient for the trivial and sign representation as f r ≡ f t − f s and f ≡ f s we cansimplify the expression further: h E ( x , k ) = 12 (cid:88) m,m (cid:48) f r | , m (cid:105) (cid:104) , m (cid:48) | + ( N − / (cid:88) l = − ( N − / (cid:88) m f l | l, m (cid:105) (cid:104) l, m | (A.7)The second term can be viewed as the electric energy of Z N as it acts trivially on thegauge 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 theHamiltonian is required. Since the calculation of these commutators for a general latticegauge theory is very lengthy, it is only sketched here. Appendix B.0.1. First order For the first order formula the ordinary commutators needto be computed. Starting with the commutator between gauge-matter interactions ondifferent links i and j , we obtain: igital quantum simulation of lattice gauge theories in three spatial dimensions H GM,i , H GM,j ]= (cid:88) x,k i λ GM ψ † m ( x ) U mn ( x, k i ) ψ n ( x + k i ) + h.c, (cid:88) y,k j λ GM ψ † m ( y ) U mn ( y, k j ) ψ n ( y + k j ) + H.c = λ GM (cid:88) x/ { boundary } ψ † n ( x + k i ) U † nm ( x, k i ) U mn (cid:48) ( x, k j ) ψ n (cid:48) ( x + k j ) − H.c. + ψ † m ( x ) U mn ( x, k i ) U † nm (cid:48) ( x + k i − k j , k j ) ψ m (cid:48) ( x + k i − k j ) − H.c. = λ GM (cid:88) x/ { boundary } U † W U W ( ψ † n ( x + k i ) ψ n ( x + k j ) − H.c. ) U † W U W + U W U † W ( ψ † n ( x ) ψ n ( x + k i − k j ) − H.c. ) U W U † W (B.1)where we used the unitary operators U W from Sec. 3.2.2 to reduce the gauge-matterterms to pure fermionic tunneling terms, thus allowing to estimate this expression: (cid:107) [ H GM,i , H GM,j ] (cid:107) ≤ λ GM N links d d U (B.2)where d U is the dimension of the representation of U under the gauge group and thereforethe operator norm of the tunneling term. In the next step, the commutator betweenthe matter- and gauge-matter interactions is calculated:[ H M , H GM,i ]= (cid:34)(cid:88) x M ( − x ψ † n (cid:48) ( x ) ψ n (cid:48) ( x ) , (cid:88) y,k i λ GM ψ † m ( y ) U mn ( y, k i ) ψ n ( y + k i ) + H.c (cid:35) =2 (cid:88) x M λ GM ( − x ψ † m ( x ) U mn ( x, k i ) ψ n ( x + k i ) − H.c. (B.3)We rewrite this expression again in terms of the unitary operators U W which allows usto bound the commutator in the following way: (cid:107) [ H M , H GM,i ] (cid:107) ≤ M λ GM N links d U (B.4)For the commutator with the electric part the whole gauge-matter Hamiltonian isconsidered as every part does not commute with H E . To bound this expression fromabove we write H GM again in terms of the unitary operators U W , similar to the previouscalculations and obtain: (cid:107) [ H GM , H E ] (cid:107) ≤ N links λ GM λ E max j | f ( j ) | d U (B.5)The last commutator is the one between the magnetic and electric Hamiltonian. Sinceevery link is contained in 2( d − 1) plaquettes, the commutator is straightforwardlyestimated as: igital quantum simulation of lattice gauge theories in three spatial dimensions (cid:107) [ H B , H E ] (cid:107) ≤ λ B λ E N links d − 1) max j | f ( j ) | d U (B.6) Appendix B.0.2. Second order For a bound on the second-order formula we need tocalculate all nested commutators. The computations of them are done in the samemanner as for the ordinary commutators, there are no additional tricks required. Sincethese calculations are very lengthy, we will just give the bounds obtained for each nestedcommutator: (cid:107) [[ H B , H E ] , H E ] (cid:107) ≤ λ E λ B max j | f ( j ) | N links d − d U (cid:107) [[ H B , H E ] , H B ] (cid:107) ≤ λ E λ B max j | f ( j ) |N l links d − d U (cid:107) [[ H E , H GM ] , H GM ] (cid:107) ≤ λ GM λ E max j | f ( j ) |N links (2(2 d − 1) + 1)4 d U (cid:107) [[ H E , H GM ] , H E ] (cid:107) ≤ λ GM λ E max j | f ( j ) | N links d U (cid:107) [[ H M , H GM ] , H GM ] (cid:107) ≤ λ GM M N links dd U (cid:107) [[ H M , H GM ] , H M ] (cid:107) ≤ λ GM M N links d U (B.7)In the last step the nested commutator among the different gauge-matter Hamiltoniansneeds to be computed: (cid:107) [[ H GM,i , H GM,j ] , H GM,l ] (cid:107) ≤ λ GM N links d d U (B.8)To obtain the error bound for the whole gauge matter interactions we need to calculatehow many times the commutator from above appears. There are 2 d different gauge-matter Hamiltonians which are implemented separately. Recalling the second orderformula, this gives rise to two partial sums over the natural numbers: d − (cid:88) k =1 (cid:107) [[ H GM,k , H GM,k +1 + .. + H GM, d ] , H GM,k +1 + .. + H GM, d ] (cid:107) + 12 (cid:107) [[ H GM,k , H GM,k +1 + .. + H GM, d ] , H GM,k ] (cid:107)≤ λ GM N links d d U d − (cid:88) x =1 x + x λ GM d U N links (2 d − (cid:18) 23 (4 d − 1) + 1 (cid:19) (B.9)Inserting all these commutators into the formulas of the total trotter error will thenresult in the bounds given in Sec. 3. References [1] Eidelman S, Hayes K, Olive K e, Aguilar-Benitez M, Amsler C, Asner D, Babu K, Barnett R,Beringer J, Burchat P et al. Physics Letters B igital quantum simulation of lattice gauge theories in three spatial dimensions [2] Peskin M E 1995 An introduction to quantum field theory (Westview press)[3] Gross D J and Wilczek F 1973 Physical Review D Physical Review D Reviews of Modern Physics Introduction to quantum fields on a lattice (Cambridge University Press)[7] Aoki S, Aoki Y, Bernard C, Blum T, Colangelo G, Della Morte M, D¨urr S, El-Khadra A X, FukayaH, Horsley R et al. The European Physical Journal C Physical review letters Reviews of Modern Physics The phases of quantum chromodynamics: from confinementto extreme environments vol 21 (Cambridge University Press)[11] Cirac J I and Zoller P 2012 Nature Physics Reviews of Modern Physics Reports on Progress in Physics Annalen der Physik Physical Review Letters Annals of physics Reviews of modern physics Physical review letters Reviews of Modern Physics Physical Review B Physical review letters Physical reviewletters Physicalreview letters Physical review letters Physical review letters Physical review letters New Journal of Physics [28] Notarnicola S, Ercolessi E, Facchi P, Marmo G, Pascazio S and Pepe F V 2015 Journal of PhysicsA: Mathematical and Theoretical Physical review letters Physicalreview letters Physical Review A Physical Review Letters Physical Review X Physical Review A Nature communications Annals of Physics Physical review letters Annals of physics Physical review letters Physical Review A Physics Letters B Physical Review A igital quantum simulation of lattice gauge theories in three spatial dimensions [43] Zohar E and Reznik B 2013 New Journal of Physics EPL (Europhysics Letters) Nuclear Physics A arXiv preprint arXiv:1510.06754 [47] Feynman R P 1982 International Journal of Theoretical Physics Quantum Information & Computation et al. Nature New journal of physics Proceedings of the American Mathematical Society Physical Review D Physical Review D Physical Review D Physical Review A Journal of Physics A: Mathematical and Theoretical Communications in Mathematical Physics Monte Carlo Methods in Quantum Problems vol 125 (Springer Science & BusinessMedia)[59] Binder K 1992 Topics in Applied Physics Lattice gauge theories and Monte Carlo simulations (World Scientific)[61] Suzuki M 1991 Journal of Mathematical Physics Journal of mathematical physics Ultracold Atoms in Optical Lattices: Simulatingquantum many-body systems (Oxford University Press)[64] Hamann S, Haycock D, Klose G, Pax P, Deutsch I and Jessen P 1998 Physical Review Letters Physical Review A R20[66] Aguado M, Brennen G, Verstraete F and Cirac J I 2008 Physical review letters101