Quantum Simulation of Gauge Theories

A general scheme is presented for simulating gauge theories, with matter fields, on a digital quantum computer. A Trotterized time-evolution operator that respects gauge symmetry is constructed, and a procedure for obtaining time-separated, gauge-invariant operators is detailed. We demonstrate the procedure on small lattices, including the simulation of a 2+1D non-Abelian gauge theory.

Evaluation of observables represented by a single, hermitian operator at a single time is straight forward. This includes most static properties. Other observables, like time-separated composite operators (e.g. parton distribution functions) are more complicated. Naively, the first nonunitary operator collapses the quantum state. One resolution introduces ancillary probe and control qubits [41], which we use here. QPE [28] has been used compute linear response [6]. Another intriguing idea uses quantum sensors to implement generating functionals [42]. A second problem arises from state contamination. In Euclidean lattice field theory, we can use imaginary-time evolution to separate states overlapping with the same operator. Real-time evolution lacks this separation, so we need novel methods to project onto desired states, perhaps through a distillation-like technique [43].
Violations of gauge invariance in the time-evolution are potentially introduced through noisy gates, digitization, and Trotterization. To remedy these, [58] proposed using oracles to project into the physical subspace. In this work, we construct a Trotterized time-evolution without gaugefixing, for arbitrary gauge theories with coupled matter fields. Gauge invariance is exactly preserved after Trotterization. We connect this procedure to the transfer matrix formalism [59][60][61][62]. The effect of noisy gates on gaugeinvariance are hardware-dependent, and we leave this for future work. The error introduced by field digitization can be treated without impacting the algorithm.
A final, overarching concern is renormalization. Work in quantum computing [63,64] and elsewhere [65,66] suggest that Trotterization can be mapped to a Euclidean action arXiv:1903.08807v1 [hep-lat] 20 Mar 2019 such that the renormalization factors may be classically computed. The transfer matrix formalism makes this a natural procedure. This paper is structured as follows: in Sec. II we describe the prerequisite registers and fundamental gates for our construction. With these, Sec. III formulates gaugeinvariant U(t) for the gauge fields. Sec. IV links this procedure to the transfer matrix formalism, and outlines how renormalization may be performed classically. Sec. V explains how nontrivial correlators may be computed. Sec. VI and VII extend the formulation to include scalar and spinor fields respectively. Sec. VIII demonstrates of our method in the non-Abelian D 4 theory and the Z 2 theory with staggered fermions. We conclude in Sec. IX.

II. PREREQUISITES
Developers of quantum algorithms assume idealized qubits, and certain unitary gates operating on those qubits (a sufficient universal gateset, for example, is HADAMARD, TOFFOLI, and CNOT). This follows a pattern in classical computing, where certain constructs (e.g., floating-point numbers and arrays) are assumed to specify an algorithm. The implication of this is that any correct implementation of the primitive constructs will do: the correctness of the quantum algorithm does not depend on whether superconducting qubits or trapped ions are used.
In this section, we describe a set of higher-level constructs that allows the simulation of a general nonabelian gauge theory. We outline how they might be implemented in practice, from the lower-level qubit/gateset primitives. Hopefully these specifications are helpful to focus the implementation of those tools that are most critical.
We will simulate a theory with a gauge group G. The central construct is the G-register, analogous to a classical variable storing an element of G. The Hilbert space of an ideal G-register is H G = CG, the complex vector space with one basis element |g for every element g ∈ G.
(Equivalently, it is the space of square-integrable functions from G to C.) The Hilbert space of multiple G-registers is constructed as tensor product, so a computer with N G-registers has Hilbert space H ⊗N G . For discrete groups, ideal G-registers are easy to implement by virtue of their finite number of elements; for instance, an ideal Z 2 -register requires one qubit. For continuous Lie groups (e.g. SU (3)), ideal G-registers cannot be made with any finite number of qubits -we must settle for approximate G-registers (typically with Hilbert spacẽ H G = span{|g } g∈G for some subsetG ∈ G) which leads to intrinsic discretization error. The tradeoff between number of qubits and the size of discretization error will be crucial to constructing efficient simulations near-term.
For the Hamiltonians we consider, a set of useful primitive gates defined on the G-register is: inversion, matrix multiplication, trace, and the quantum Fourier transform. The inclusion of matter fields requires additional multiplication and Fourier transform gates, and an inner product gate. A sample implementation for these operations is specified in Appendix A for the D 4 -register case.
The inversion gate acts on a single register to produce the register of its inverse group element. This is defined in the fiducial basis by (1) The other group operation is matrix multiplication, which requires a two-register gate, defined by Here we have U × implementing left multiplication, in the sense that the content of the target register was multiplied on the left. This is all that is required here. In general, the combination of U −1 and leftacting U × permits right multiplication: Most quantum gatesets -certainly all gatesets proposed for practical use -contain the inverse of every gate as a primitive. This makes it easy, given a circuit implementing one unitary, to obtain a circuit implementing its inverse. Therefore, if U † × is available, U −1 is readily obtained with the assistance of an ancillary register.
The trace of a plaquette appears in our gauge Hamiltonian, and so to perform this operation we combine with the matrix multiplication gate with a single-register trace gate: In our construction, the final gate required on the Gregister is the Fourier transform gate [67][68][69] U F . This gate acts on a G-register to rotate into Fourier space. It is defined by The second sum is taken over ρ, the representations of G, andf denotes the Fourier transform of f . This gate diagonalizes what will be the 'kinetic' part of the Trotterized time-evolution operator. After application of the gate, the register is no longer a G-register but aĜ-register. The G-register suffices for the representation of the gauge field. In order to represent scalar matter fields, we introduce the C N -register, where N is the dimension of the fundamental representation 1 of G. This register stores a component of C N , and we will define gates allowing transformation under elements of G. As with the Gregister, it is only possible to implement a 'truncated' C N register. (For fermionic matter fields, the C N -register is unneeded, and no truncation is necessary.) Three gates are needed to make the C N -register useful in a quantum simulation. The fundamental (or adjoint) multiplication gate U ×,f is defined by The inner product gate, necessary for the field-theory kinetic term (sometimes called the 'hopping' term, to distinguish it from the quantum mechanical kinetic term), is defined by The inner product is generally a complex number, but this gate extracts only the real part while remaining unitary.
Happily, this is all required for simulation of field theory, where only the linear combination φ † 2 φ 1 + φ † 1 φ 2 appear. Lastly, the scalar Fourier transform gate is defined analogously to U F above, except for C N : we notate this U C N F . The resulting register is a C N -register. Since the quantum Fourier transform on C N factorizes into N separate transforms on C, the implementation of U C N F consists of N separate calls to U C F .

III. PURE GAUGE
In this section, we construct a time-evolution operator suitable for the simulation of pure gauge theories on a quantum computer. We take care to exactly preserve gauge invariance (up to unavoidable gate errors and noise). We take the standard approach, by discussing the Hamiltonian formulation of the theory, and then Trotterizing the Hamiltonian into a sum of pieces easily implemented via the primitive constructs of Sec. II. This approach will be extended in Sec. VI and VII to include scalar and fermionic matter fields, respectively.
We begin by reviewing the Hamiltonian formulation of lattice gauge theory [70]. Each gauge link has a local Hilbert space CG. Combined, the L links on our lattice have Hilbert space H = CG ⊗L .
On a lattice with N sites, the space of gauge transformations is G N , where one element of G must be specified at each site. The transformation rule for a link U ij from Gauge invariance demands that two states differing only by a gauge orbit be considered physically equivalent. The physical Hilbert space 2 , then, is H P = CG ⊗L /φ(G ⊗N ). 2 The Gribov ambiguity prevents us from unambiguously assigning a single representative configuration to each gauge orbit; however, we do not take that approach. Physical states in the Hilbert space are not configurations, but formal linear combinations of configurations, and these can be constructed unambiguously.
H P can be obtained from the larger space H by a gauge-symmetrization operator, which acts as a projection operator: This discussion of the gauge-invariance of physical states is equivalent to the imposition of Gauss's law. An implementation of the operator P for U (1) is given by [58]; however our method does not require implementing P .
We define a Hamiltonian which acts on the entire space H of the gauge theory, implicitly defining a 'physical' Hamiltonian on the subspace H P . As will be made clear in later sections, this choice of Hamiltonian is motivated by a desire to recover the Euclidean Wilson action.
The first term, H V , is a sum over spatial plaquettes p and can be thought of as a potential term. The trace is taken in a suitable representation of the gauge group G -for matrix groups, this is typically the fundamental representation. The second term, H K , is the quantummechanical kinetic term describing the mechanics of a free particle moving on the surface of G. Therefore, π 2 ij is the Laplace-Beltrami operator on the surface G, and π in is the momentum operator conjugate to U ij . This prescription gives us no guidance for how to form the kinetic term of a discrete gauge theory -this is addressed via the transfer matrix formalism in Section IV.
Not only is the Hamiltonian gauge-invariant, but each term is gauge-invariant. This is crucial, as it allows us to Trotterize without breaking gauge invariance at any order in ∆t. A time evolution operator suitable for implementation on a quantum computer is then Below, we will use the notation U(∆t) to indicate a single ∆t step. We now describe how the operators e −iH K ∆t and e −iH V ∆t may be obtained. First, each part (kinetic and potential) of the Hamiltonian consists of mutually commuting, local terms, and each such term may be treated individually and sequentially without approximation. That is, The plaquette evolution U (1) V is constructed by first preparing the product of the plaquette in an ancillary register, and then operating on the ancillary with U Tr ( 1 g 2 a ). U K is implemented by diagonalizing via U F , and then applying a diagonal unitary: Fig. (1). The implementation of U phase depends strongly on the group and the implementation of the G-register, but in all cases it is a diagonal operator. The total gate requirements for the propagation are shown in Table I.
When preparing an initial state, it is important to ensure that it lies within the physical subspace H P . One procedure for projecting onto the gauge-invariant states for an Abelian gauge theory is provided in [58].
Although the preparation of the ground state |Ω is beyond the scope of this paper, it is worth discussing that common proposals typically require the time-evolution operator U(t) as constructed in this section. Assuming such methods are used, nonequilibrium physics can be computed by preparing the ground state with U(t), and then evolving with a different time-evolution operatorŨ(t) [11]. A typical expectation value of interest, accessible in such a scheme, is whereÊ measures the expectation of a plaquette (i.e. the energy density of the gauge field). Other expectation values of interest are discussed in Sec. V below.

IV. THE LATTICE PATH INTEGRAL
Substational insight into QCD has been achieved via Euclidean-spacetime lattice methods on classical computers. Accordingly, the Euclidean lattice path integral occupies a dominant place in the QCD literature and the minds of practitioners. This is in contrast to quantum simulations where the Hamiltonian formulation appears more natural. In this section, we connect the Hamiltonian formalism to the Euclidean lattice path integral via a transfer matrix. This implies an important practical consequence: the determination of appropriate bare coupling constants may be performed on a classical computer. This procedure also guides us in constructing a Hamiltonian for a discrete gauge theory (discussed in [71]).
Our formalism for the transfer matrix, T , slightly differs from the usual one [59][60][61][62]72] in that T is defined on the entire space H. The usual transfer matrix is defined only on the physical Hilbert space, and may be obtained by projection with P . In the course of performing this projection, the lattice path integral with the Wilson action is recovered.
The transfer matrix is an approximation to imaginarytime evolution by one unit of time: T ≈ e −H . We construct the transfer matrix in the fiducial basis of H from separate kinetic T K and potential T V contributions via the split-operator approximation where the potential term resembles the product of spatial plaquettes as appears in the Wilson action (15) Here a is the spatial lattice spacing, and a 0 the temporal spacing; we do not assume an isotropic lattice. We have borrowed from lattice field theory the Wilson plaquette and µ, ν are restricted to space-like directions.
For the kinetic term there is no coupling between different links, but it is related in lattice field theory to plaquettes containing timelike links: Note that, like the Trotterized real-time evolution, [T, P ] = 0 due to the fact that T V and T K individually commute with P . When the transfer matrix is restricted to physical states, it exactly reproduces the usual lattice path integral, with the Wilson action.
Here Tr P denotes the trace over only the physical subspace H P . The inverse temperature is denoted β, and Z is the Euclidean lattice approximation to the partition function. Now that we have a transfer matrix corresponding to the Wilson action, we may connect it to the Hamiltonian. This is done in detail by Creutz [72] for the gauge-fixed transfer matrix, so we will only summarize the gaugefree result (which is essentially the same) here. The potential part T V of the transfer matrix is exactly e −H V . The relationship between the kinetic part H K of the Circuits for the propagation of a pure-gauge lattice field theory. The first circuit implements U K on four links (in general, L links are needed). The second circuit shows the application of U (1) V to a single plaquette Re Tr U † 13 U † 34 U24U12, and must be applied to every plaquette in the theory. Note that in these circuits, we use a doubled line to represent a G-register, rather than a single qubit.
Hamiltonian and T K is more subtle. It is shown in [72] that, in the limit where the temporal lattice spacing a 0 is taken to 0, the kinetic part of the Hamiltonian is recovered.
With a mapping between the Hamiltonian and lattice action through the transfer matrix, the bare constants in Eq. (9) for a desired renormalized, physical set of couplings are the same as the lattice action ones. This is particularly important for ensuring that Lorentz invariance is recovered in the continuum limit, but also allows us to perform simulations at known values of physical parameters. In particular, in a lattice QCD simulation, the task of determining appropriate bare parameters for a desired physical pion mass may be performed cheaply on a classical computer, reserving the quantum processor for the difficult real-time evolution.
Clearly, the Hamiltonian evolution on a quantum computer is closely related to standard classical lattice simulations. Indeed, the two differ only in: the discretization of the gauge group (quantum computers have smaller G-registers), whether the Trotterization is done in real or imaginary time, and the exact form of the kinetic term. The first difference can be removed by performing a classical simulations with reduced precision, and the others are reduced at small temporal lattice spacing.
In the case of the Trotterization ∆t being equal to a, the quantum evolution is simply an analytic continuation of the usual isotropic lattice calculations performed in lattice QCD. This implies that Lorentz invariance is recovered in the continuum limit.
We conclude by exploiting this procedure to produce a Hamiltonian for a discrete gauge theory -that is, a gauge theory where the gauge group is a discrete group (such as D 4 , demonstrated below), rather than a connected manifold. In this case, the kinetic part of the Hamiltonian is determined not by a Laplacian on a manifold, but can be derived from relevant portion of the transfer matrix by taking a logarithm.
where T V and T K are defined exactly as in Eq. (14). Note that H K is only well-defined if T K is positive-definite: this fact is established for a general discrete group (with the Wilson action) in Appendix B. As a consequence of Schur's lemma, the kinetic transfer gate defined here is diagonal in the Fourier basis. This implies that, where H K is obtained in this fashion, U K may be easily implemented from U F .

V. CORRELATORS
Having constructed time-evolution for a pure gauge theory, we now must show how observables of physical interest may be obtained. In this section, we show how to obtain: a plaquette-plaquette correlator, and the expectation value of an arbitrary, temporally extended Wilson loop. We accomplish both by the technique set out in [41]. To begin, we recall how one can compute an expectation value of a unitary operator U [73] in a given state |Ψ . Introducing a single ancillary qubit, we construct a controlled unitary operator U C , defined by U C |Ψ |0 = |Ψ |0 and U C |Ψ |1 = U |Ψ |1 . (19) (Standard quantum gatesets make the construction of this operator straightforward.) Generally, the expectation value of U has both real and imaginary parts. These must be measured separately. A measurement of Re Ψ| U |Ψ is given by initializing the ancillary qubit to |+ ≡ 1 √ 2 (|0 + |1 ), performing evolution via U C , and then measuring σ x on the ancillary qubit.
For Im Ψ| U |Ψ one performs the same process except ending with a measurement of σ y .
With this procedure in mind, we can show how to compute a correlator of the form In QCD this Wilson plaquette correlator corresponds to a glueball propagator. The operator in Eq. (21) is clearly not unitary, so it cannot be directly evaluated by means of the procedure described above. Instead, using the Hermiticity of W µν (x), we construct a parameterized family of unitary operators. The derivatives of this family give the operator of interest. Explicitly, this is done via a time-dependent perturbation of the Hamiltonian: Time evolving forward in time with H 1, 2 , and back with H 0 , allows us to measure the expectation value C( 1 , 2 ) ≡ Ψ| U(−t)U 1, 2 (t) |Ψ . Differentiating twice with respect to the perturbation yields the plaquetteplaquette correlator as desired.
In practice, this derivative must be obtained by finite differencing, after computing C( 1 , 2 ) for several small values of . This approach is naturally extended to the computation of k-time correlation functions, at the cost of requiring the kth numerical derivative.
The terms W µν (x) added in the perturbation also appear in the original Hamiltonian, and are readily implemented with the primitive gates defined in Section II.
We can use the same method to compute the expectation values of temporally-extended Wilson loops, where a, b are color indices. We can use this to derive a perturbed Hamiltonian which can be utilized in a correlator C ab ( 1 ,˜ 1 , 2 ,˜ 2 ) ≡ Ψ| U(−t)U ab (t) |Ψ We now take the difference of two second derivatives of the resulting correlator: Summing over the a, b yields the desired, gauge-invariant trace. The imaginary part can be obtained similarly. For many gauge groups, the different correlators (U † ij (t)) ba (U ij (0)) ab are related by gauge symmetry. Consequently, for these gauge groups, the correlator in Eq. (25) need only be evaluated for one particular selection of a and b, as all the other terms will be equal. In particular, this is true for D 4 (simulated below) and SU (N ) in the fundamental representation.
The perturbation introduced, unlike the one used for plaquette-plaquette correlators above, is not gaugeinvariant, meaning that the state during this time evolution does not lie within the physical subspace H P . Also unlike the plaquette-plaquette perturbation, this perturbation contains terms not present in the original Hamiltonian, and requires gates not specified in Section II. The gates required are diagonal (in the fiducial basis) phase gates on individual G-registers.
This method can also compute a non-gauge-invariant expectation value, e.g. Tr U ik (t)U ij (0) where j = k. Provided the initial state is correctly prepared, this expectation value must be 0.

VI. SCALAR FIELDS
Here we extend our Hamiltonian construction to include scalar matter fields which transform in a D-dimensional representation of G. The Hilbert space defined in Section III is extended via tensor product to include D scalar degrees of freedom at each site. The new Hilbert space is H = (H G ) ⊗L ⊗(H S ) ⊗N , where H G is the one-link Hilbert space defined above, N is the number of sites, and H S is the one-site Hilbert space. For H S , it is convenient to work in the basis of occupation number, writing H S = span{|n 1 · · · n D : n i = 0, 1, . . .}.
As before, we first consider how gauge transformations act on H. A gauge transformation is still specified by The action of a gauge transformation V on the Hilbert space is (26) The gauge-symmetrization operator is still defined by Eq. (8), but acts on the expanded Hilbert space.
To represent the physical Hilbert space on a quantum computer, we now make use of the C N -registers transforming under G. For the scalar fields we consider, the Hamiltonian is Although we assume scalar interactions only to order φ 2 , the generalization to higher-order couplings is trivial and introduces no complications into what follows. We construct a Trotterized time evolution operator, just as in Section III, by decomposing into potential and kinetic terms. The new gates required for the scalar fields are U  Uij φi) , and (29) Note that the mass and hopping gates commute, as both are diagonal in the fiducial basis. The same-site φ † i φ i terms generated by the hopping term of the Hamiltonian have been absorbed into the mass gate, assuming a square lattice of dimension d. These gates, together with the pure-gauge U (1) V,K , are combined to form With these, the full time-evolution operator for a single time-step is A single time step, therefore, consists of five steps on the quantum computer: 1. For each plaquette, compute the plaquette product and apply U Tr (∆t/g 2 a). This implements U V,gauge .
2. For each link, apply the scalar hopping gate U K,scalar . This is the general algorithm for time-evolving gauge fields in the presence of matter. Before moving on to discuss scalar field correlators, we note that the transfer matrix discussion of Sec. IV could be repeated here, including the scalar fields, to connect the evolution of the quantum computer to a classical, Euclidean spacetime path integral. Thus, renormalization calculations again may be performed classically.
Steps 1 and 4 above are the same as for a pure-gauge simulation.
Step 2 requires two calls to U ×,f , two calls to U −1 , and one call to U (·,·) . Step 3 is a registerimplementation dependent diagonal phase gate.
Step 5, analogous to the implementation of U (1) K,gauge , entails 2N calls to U C F and N diagonal phase gates per site, where N is the dimension of the scalar's representation.
The procedure for computing plaquette correlators is the same as in Section V.
The construction of Section V may be generalized to obtain the gauge-invariant two-point function of a fermion, connected by a Wilson line: ψ (y, t)W ψ(x, 0) . However, for a general Wilson line, this construction requires one derivative to be taken at every timeslice. As the derivatives must be taken via finite differencing, this is not a plausible method for computing these correlators, even in the far future. Constructing a more efficient method remains an open problem.

VII. FERMION FIELDS
Here, we consider the inclusion of spinor matters fields, which can be done precisely the same as the the method described above for scalar matter. In the implementation of spinor fields on a quantum computer, while the register is trivial to specify and exact, a significant technical difference emerges: the fundamental operators provided on most digital quantum computers are bosonic, meaning that operators acting on different qubits will commute with each other. For spinor fields, we want fermionic operators, which anticommute at different physical sites. We will translate commuting operators into anticommuting operators via the Jordan-Wigner transformation, paying the price that local fermionic operators are generated by nonlocal bosonic operators. More efficient, albeit complex, schemes, relating local fermionic operators to more local bosonic operators have been devised [14][15][16].
The Hilbert space of a gauge theory with d-component fermions can be written where there are L links, N sites, and H F = span{|0 , |1 } is the Hilbert space of a single fermion on a site. Note that, as a consequence of the exclusion principle, H F can be put on the quantum computer without approximation or truncation.
Writing the Hamiltonian using operators that obey the standard fermionic anticommutation relations is easily accomplished in the absence of gauge fields via the Jordan-Wigner transformations [13]: each fermionic operator pair is given an index from 1 · · · dN , and the ith operators are constructed from bosonic operators a i and a † i : In the presence of gauge fields, we must ensure that the fermionic operators obey the appropriate gauge transformation law. To be concrete, consider when the fermions are transform in the fundamental of SU (3). Then at each site i, for each spinor index α, there are three annihilation and creation operators, which must transform under gauge transformations as where a, b are color indices. There are now three times as many fermion operators, as for each site and each spinor, there must be three different fermionic degrees of freedom. This may be accomplished if the bosonic operators a, a † from which the fermions are constructed transform in the fundamental of SU (3). Thus each operator has a color index, and a a → g ab a b . Crucially, no change to the Jordan-Wigner procedure is need since the σ-matrices are independent of the gauge group (however the σ matrices do inherit the color index).
With fermionic operators constructed, the timeevolution operator is constructed in a similar fashion to that from Section VI above.
The key difference is that the spinor kinetic and potential operators are derived from the Dirac equation. The governing Hamiltonian is now and the single-site and single-link evolution operators are defined as before: K,spinor (i, j) = e −i∆tψj Uij ψi (38) An implementation of Z 2 fermions is demonstrated in Section VIII.

VIII. DEMONSTRATION
In this section we present two concrete demonstrations of the method described above: first, a two-plaquette theory with the discrete, non-Abelian gauge group D 4 , and second, a Z 2 gauge field coupled to fermions. For all our simulations, we work in units where the lattice spacing is a = 1. We perform classical simulations of a quantum computer using qiskit [? ? ], without simulated noise in order to demonstrate the correctness of the implemenation.
We simulate the D 4 gauge field on the two-dimensional lattice shown in Fig. (2). This simulation requires a four D 4 registers, and uses a total of 14 qubits: 12 for physical degrees of freedom, and 2 ancillary qubits. Note that, for brevity, we have broken with the notation of previous sections, in referring to a link not by the source and sink sites, but instead with a single direct index 0 . . . 3.
We define a trace on D 4 (not uniquely specified by the group structure) by embedding D 4 into U (2), and defining the trace via the fundamental representation of that Lie group. The embedding of D 4 < U (2) is generated by the elements σ x and iσ z . The lattice action for this model is The last term is a sum over all temporal plaquettes on the 2 + 1 lattice. The resulting Hamiltonian is where log T (1) K is the one-link kinetic term, determined as discussed in Section IV.
The starting state of our D 4 simulation is the gaugeinvariant projection of the eigenstate of the operators U i , where all links have been assigned to the identity matrix. We evolve the system for t = 10 with a Trotterization step of ∆t = 0.2. The circuits used are detailed in Appendix A. In total, the quantum simulation of Fig. (3) entailed ∼ 200 entangling gates per Trotterization time step. This is roughly the resources recently used to compute the ground state energy of the water molecule [74]. Fig. (3) shows the average plaquette energy of the left plaquette, as a function of time. The only systematic error is Trotterization, as there is no noise on the simulated quantum computer.
For the same lattice gauge theory, we also demonstrate the procedure described in Section V for determining the expectation value of a Wilson loop with temporal extent t: Re Tr U † 0 (t)U 0 (0) . Simulated results are shown in Fig. (??). For the finite differencing, we perform simulations with ( 1 , 2 ) = (0.1, 0.0), (0.0, 0.1), (0.1, 0.1), and the same for the parameters˜ 1 ,˜ 2 . As before, a Trotterization time step of ∆t = 0.2 is used. We note that there is a tradeoff between the size of (smaller deformations improve the finite differencing calculation) and the number of samples required. In order to test our technique with matter fields, we simulated staggered fermions transforming in a Z 2 gauge field, in 1 + 1 dimensional spacetime. A similar model was also simulated on a quantum computer in [10]. The lattice Hamiltonian for this model is where the indices i, j denote lattice sites, the sum is taken over all pairs of adjacent sites, and χ is a one-component fermionic operator.
Using the Jordan-Wigner procedure, we translate from the fermionic operators χ i into a set of bosonic operators σ i : for site i. In one-dimension without periodic boundary conditions, this results in a local bosonic Hamiltonian.
There are two distinct families of bosonic operators: the σ i which live at a site and correspond to fermions, and the σ ij on the link from site i to site j, corresponding to the Z 2 gauge field. The Trotterization of this Hamiltonian takes four steps, corresponding exactly to the four terms in Eq. (42). In this way full translational invariance (neglecting the boundary) is maintained even in the approximated evolution. The starting state of this Z 2 simulation is obtained by projecting the |0000000 state to the gauge-invariant space, and then applying the gauge-invariant excitation operator χ † 1 χ † 2 σ (z) 12 . The number density of the left-most site as a function of t is shown in Fig. (5) where the simulated results are compared to the exact values. Since we do not begin in an eigenstate of the Hamiltonian, this is not constant. The difference between the two is due to the Trotterization, and as with the D 4 simulation may be made systematically smaller by reducing the lattice spacing. This simulation uses 8 qubits: 4 fermion sites, 3 gauge links, and a single ancillary. For this lattice, 36 CNOT gates are used per timestep, putting the simulation (slightly) out of reach of present-day processors. A timestep of ∆t = 0.6 is used for T = 30.

IX. DISCUSSION
In this work, we have given a general procedure to construct a unitary time-evolution operator and timelike separated correlation functions for use on an idealized universal quantum computer. This construction ensures that the Trotterization never mixes the physical and unphysical states within the Hilbert space, exactly preserving gauge invariance and alleviating some of the concerns expressed in the literature. Further, we show how gaugeinvariant correlators (including a temporal Wilson loop) are constructed. We have demonstrated the efficacy of the procedure by the implementation of two gauge theories on a quantum simulator and a quantum computer. Additionally, we have discussed the classical procedure for determining the bare couplings to be used in a quantum simulation.
Critically, we have defined a set of primitive constructs which may be used in the construction of gauge theory simulation. The practical and efficient implementation of such registers and gates is an important open question in the literature. Finally, future work must address the issue of how the behavior of different noisy gates affects the mixing of the physical and unphysical state, whether this can be mitigated by special constructions or errormitigating codes within our construction, and indeed whether such mitigation is necessary. The first matrix represents a π 2 rotation, and the second represents reflection.
The D 4 register is implemented with three qubits. The state |abc corresponds to the matrix (A2) Thus, the two least-significant bits specify an amount of rotation, to be followed by a flip of the most-significant bit is 1. This establishes the isomorphism needed between CD 4 and the three-qubit Hilbert space on a quantum computer. It remains to describe the inversion, multiplication, trace, and Fourier transform circuits. The inversion FIG. 6. Classical-inspired quantum circuits for D4 inversion (top) and multiplication (bottom). The least-significant qubits are shown on the top. For the multiplication circuit, the target register is in the bottom three qubits.
and multiplication circuits are both equivalent to classical circuits (as they consist exclusively of controllednot gates), and are shown in Fig. (6). The trace circuit is a three-qubit controlled phase gate, defined by U Tr (θ) |000 = e −2iθ and U Tr (θ) |010 = e 2iθ (with all other states picking up no phase). Finally, the Quantum Fourier Transform which diagonalizes the D 4 momentum operator is given by (A3) Rather than attempt to determine a quantum circuit by hand (as general methods exist [75][76][77]), we perform a classical simulated annealing search to find a short circuit that implements exactly this unitary. The circuit found is shown in Fig. (7). The annealing algorithm will be detailed in a future work. Since the size of F does not scale with volume, this method does not suffer from an exponentially large physical Hilbert space. After being diagonalized by the Fourier transform, the kinetic gate corresponds simply to a phase gate on the most-significant qubit, along with a phase gate on the state |000 :