Quantum simulation of gauge theory via orbifold lattice

We propose a new framework for simulating U(k) Yang-Mills theory on a universal quantum computer. This construction uses the orbifold lattice formulation proposed by Kaplan, Katz, and Unsal, who originally applied it to supersymmetric gauge theories. Our proposed approach yields a novel perspective on quantum simulation of quantum field theories, carrying certain advantages over the usual Kogut-Susskind formulation. We discuss the application of our constructions to computing static properties and real-time dynamics of Yang-Mills theories, from glueball measurements to AdS/CFT, making use of a variety of quantum information techniques including qubitization, quantum signal processing, Jordan-Lee-Preskill bounds, and shadow tomography. The generalizations to certain supersymmetric Yang-Mills theories appear to be straightforward, providing a path towards the quantum simulation of quantum gravity via holographic duality.


Introduction
Quantum simulation is among the most important applications of quantum hardware, both for near-term and fault-tolerant quantum computation, addressing the capability of quantum devices to probe, calculate, and simulate real problems appearing in the physical world.Theoretically, constructing efficient quantum simulation algorithms enhances support for the claim of the quantum Church-Turing Thesis, which states that one can simulate all physical processes using quantum circuits with reasonable costs in time.Practically, quantum simulation of fundamental physical processes may shed light on complex phenomena appearing in quantum gravity, quantum cosmology, subatomic particle physics, condensed-matter physics, cold-atomic physics, and statistical physics where classical computers encounter insurmountable challenges.We are in an era of quickly developing quantum technology, where near-term quantum computers may perform tasks surpassing the current capabilities of classical computation.In the long-term, it is reasonable to expect that quantum devices will perform universal faulttolerant quantum computation, allowing us to run quantum algorithms reliably (see some general introduction about this topic in [1][2][3]).
Out of the many applications of quantum simulation, implementing quantum field theory is a particularly promising area for realizing a useful quantum advantage.Already, detailed studies have been carried out demonstrating the efficiency of simulating the dynamics of scalar field theories [2,[4][5][6].Many outstanding problems in quantum field theory addressable by quantum computation concern the properties of lattice gauge theories.Gauge theories are not only the foundation of particle physics and condensed-matter physics, but also play a critical role in quantum error correcting codes and topological quantum computation.For instance, the celebrated toric code developed by Kitaev is naturally understood in terms of a Z 2 lattice gauge theory.Much attention has been given recently to the task of simulating lattice gauge theory on a quantum computer.Most proposals use the traditional Hamiltonian formulation of lattice gauge theory developed by Kogut and Susskind (KS) [11].Yet another Hamiltonian formulation is readily obtained from the lattice Yang-Mills action based on the orbifold construction introduced by Kaplan, Katz, and Unsal (KKU) in 2002 [12].In this paper, we will study how the orbifold construction can be used for digital quantum simulation and analyze its advantages and disadvantages compared to an approach based on the KS formulation.
The original motivation for the orbifold construction was to realize supersymmetry on a lattice; by performing an orbifold projection on a supersymmetric Yang-Mills Another promising approach is to start directly from a finite-dimensional model which possesses a continuum limit in the universality class of the desired quantum field theory [7][8][9][10].matrix model, one obtains a lattice theory preserving a few supersymmetries, which is (for (1 + 1) and (1 + 2) dimensions) sufficient for the complete restoration of all supercharges in the continuum limit.This idea was subsequently applied to Euclidean theories [13][14][15] and inspired alternative approaches to realizing supersymmetry on a Euclidean lattice, e.g.Refs.[16][17][18][19][20].In these formulations, no parameter fine tuning is needed to achieve the supersymmetric continuum limit to all orders in perturbation theory.Numerical results [21][22][23][24] further support the expectation that fine tuning is not necessary even at the nonperturbative level.
Therefore, the orbifold construction has a clear advantage when applied to supersymmetric theories.As far as Euclidean theories are concerned, there may be no clear motivation for applying the orbifold construction to non-supersymmetric theories, as traditional lattice regularizations such as Wilson's plaquette action [25] are well understood and already sufficiently useful.However, in this paper, we will show that the situation is different for real-time quantum simulation.In fact, the orbifold construction could provide important, non-substitutable inputs when simulating lattice gauge theories in a quantum computer.Although we consider U(k) gauge theory as a specific example, we expect that our method can be generalized to other gauge groups useful for, for instance, studies of the Standard Model in particle physics.
The original motivation of the orbifold lattice construction was the study of supersymmetric Yang-Mills theory, and more specifically, quantum gravity via holographic duality.It appears that the construction we discuss in this paper can be generalized to supersymmetric theory.Therefore, together with a recent paper [28], which studied the approach based on the matrix models, this paper may serve as the first step towards quantum simulation of deep problems in quantum gravity via holography, such as the black hole information puzzle and emergent spacetime.
This paper is organized as follows.In Section 2, we review the orbifold construction and discuss its Hamiltonian version.In Section 3, we show how to realize this model as quantum simulation, particularly for the task of preparing ground states and measuring observables.In Section 4, we compare the orbifold and Kogut-Susskind approaches as platforms for quantum simulation.In Section 5, we present some topics for future work.In Appendix A, we explicitly show how the lattice theory we consider is obtained from an orbifold projection on a matrix model.In Appendix B, we review the Kogut-Susskind formulation of lattice gauge theory.Appendix C contains some notes on an alternative digitization scheme, which may be more practical than that considered in the main text in some cases.

Orbifold construction of lattice gauge theory
In this section, we introduce the orbifold construction of pure Yang-Mills theory on a lattice.The field content and Lagrangian are explained in Sec.2.1.The orbifold construction uses non-compact variables rather than compact variables (unitary link variables); in Sec.2.2, the relation to the formulation with unitary link variables is made clear.The Hamiltonian formulation, which is used for the implementation on a quantum computer in later sections, is introduced in Sec.2.3.The symmetry of the orbifold lattice at the discretized level is examined in Sec.2.4.
The adjective 'orbifold' comes from the original construction [12], which obtained the lattice action from a matrix model via the orbifold projection.We review the details of this construction in Appendix A.
For concreteness we consider the (3 + 1)-dimensional theory.Essentially the same construction works for (2 + 1)-and (1 + 1)-dimensional theories as well, as we will briefly see in the end of Sec.2.1.

Orbifold lattice
The 'orbifold lattice' version of U(k) Yang-Mills is given as follows.We introduce a label of lattice points ⃗ n = (n x , n y , n z ), where n x , n y , n z = 1, 2, • • • , L. Then the 'orbifold lattice' is the gauged matrix quantum mechanics with the gauge group ⃗ n U(k) ⃗ n which contains k × k complex matrices x ⃗ n , y ⃗ n and z ⃗ n living on the links connecting ⃗ n and ⃗ n + x, ⃗ n + ŷ and ⃗ n + ẑ, respectively.As we will see shortly, the unitary link variables come out of these complex matrices.We use a bar to denote Hermitian conjugate, i.e., x = x † .Then x⃗ n , ȳ⃗ n and z⃗ n are regarded as the link field with the opposite direction, i.e., from ⃗ n + x, ⃗ n + ŷ and ⃗ n + ẑ to ⃗ n, respectively.We also add the gauge field A ⃗ n living on each site ⃗ n, which will be identified with the temporal component of the gauge field The orbifold construction played a key role to find a supersymmetric lattice.For the purpose of this paper, the use of the orbifold projection is not crucial; just by accepting the action given in Sec.2.1 as a starting point, it is possible to understand how the orbifold lattice serves as a lattice regularization of Yang-Mills theory.
In group theory language, x ⃗ n , y ⃗ n and z ⃗ n are bi-fundamental representations of U(k in (1 + 3)-d theory.The Lagrangian is given by (2.1) Here we have used the notation |M | 2 = M M † for any matrix M .The trace Tr is over a k × k matrix.The covariant derivative D t is defined by A local U(k) gauge transformation at site ⃗ n parametrized by Ω ⃗ n is given by and If we take the A = 0 gauge, then the Gauss-law constraint coming from ∂L ∂A = 0 is given by (2.5) Here µ = 1, 2 and 3 stand for x, y and z, respectively.Note that the coupling constant is written as g 1d because this is the same as the coupling constant of the matrix model (which can be regarded as '1d QFT'), which is used to obtain this action via the orbifold projection; see Appendix A for details.The relation between this coupling and the coupling of the 4d theory will be shown shortly.
Though we have introduced 'lattice points' ⃗ n, it is not a legitimate lattice field theory yet; even the 'lattice spacing' is missing, and there is no kinetic term!Rather, at this moment, it should be called a quiver matrix model.To make a quiver matrix model into a lattice field theory, we use dimensional deconstruction [29].The basic idea of the dimensional deconstruction is that, when some fields in a quiver theory have nontrivial vacuum expectation values, a lattice structure for extra spatial dimensions can be generated effectively.In the case of the orbifold lattice, we assume that the fields x, y and z are fluctuating around 1 √ 2ag 1d • 1 k , allowing us to identify a lattice spacing a (we will see how to justify this assumption shortly).Under this limit, writing x, y, z and A t as where s 1,2,3 and A 1,2,3 are Hermitian, and we obtain up to O(a) corrections.In this way, we can obtain a lattice regularization of (3 + 1)-d YM theory coupled to three scalar fields s 1,2,3 .Actually, the form (2.6) and the requirement that x, y, z ≃ • 1 k are too restrictive.More precisely, we only need to have to ensure the right continuum limit.Here, ≃ means the equality up to terms negligible as a → 0. Note that such small fluctuations correspond to the scalar fields, as we will see shortly.
Note that the condition (2.9) is gauge-invariant, while x, y, z ≃ 1 • 1 k is not.The following parametrization [30] gives us a convenient way to obtain the continuum The coupling constant g 2 4d in (2.7) is the bare coupling defined at the cutoff scale.
theory in a gauge-invariant manner: 2ag 1d e a 5/2 g 1d s 1 e ia 5/2 g 1d A 1 , 2ag 1d e a 5/2 g 1d s 3 e ia 5/2 g 1d A 3 . (2.10) Clearly, the condition (2.9) is violated when the scalars s 1 , s 2 and s 3 deviate too far from zero.This is known as the moduli stabilization problem [12].Modulo this problem, the orbifold construction gives an alternative to more traditional lattice regularizations.
We would now like to apply this construction to pure YM theory without the scalars s 1 , s 2 and s 3 .To accomplish this, we simply add (2.11) to the action in (2.1).In the continuum, this is nothing but the scalar mass term: With such a mass term the expansion about xx ≃ y ȳ ≃ z z ≃ 1 2(ag 1d ) 2 • 1 k is justified, resolving the moduli stabilization problem.Furthermore, by taking m 2 large enough, we can completely eliminate the scalars.In the original references of the orbifold construction, this mass term could not be used because the main motivation was a supersymmetric lattice theory, and the scalar mass breaks supersymmetry.In our case, nothing forbids us from adding this scalar mass term since we are not interested in either supersymmetry or the scalar fields, for now.If, however, we consider the simulation of supersymmetric theories, then the moduli stabilization problem comes back.In this case there are other ways to handle it; see Refs.[22,[31][32][33] for resolutions.
The biggest difference from the Wilson's plaquette action is the use of the noncompact variables.For the orbifold construction, the gauge-invariant path-integral measure is the flat measure on R 2k 2 , where µ,⃗ n and x µ,⃗ n are the real and imaginary parts of x µ,⃗ n .On the other hand, in the Wilson's plaquette action, the link variables are unitary, and the integral is defined using the Haar measure on the group manifold U(k), which is compact.Note also that, because of (2.7), in the continuum limit (a → 0, g −2 4d ∼ − log a), the coupling constant g 2 1d scales as − 1 a 3 log a .Therefore, the continuum limit corresponds to the strong-coupling limit in the original matrix model.
In the orbifold construction, the U(1) part always appears by construction.Hence it is impossible, or at least there is no known way, to construct the SU(k) orbifold lattice.However, it is not necessarily a problem.If all the fields are in the adjoint representation such as pure Yang-Mills or maximal super Yang-Mills, the U(1) part is free and decoupled from SU(k), so the local dynamics in the SU(k) sector is not affected at all.It is not the case e.g., when matters in the fundamental representation are introduced, but still, the U(1) sector can decouple from the low-energy dynamics if the U(1) is not asymptotically free.
The generalizations to other dimensions are straightforward.For the (2 + 1)dimensional theory, we use ⃗ n = (n x , n y ) to label the lattice points and the link variables x ⃗ n and y ⃗ n .The lattice Lagrangian is obtained by omitting z and z from (2.1).The coupling constant is g 2 3d = a 2 g 2 1d .For the (1 + 1)-dimensional theory, we use n = n x as the label for the lattice points, and the link variable x n .The lattice Lagrangian is obtained by omitting y, ȳ, z and z from (2.1).The coupling constant is g 2 2d = ag 2 1d .

Connection to the unitary-link formalism
Unlike the orbifold construction, the Kogut-Susskind formulation uses unitary link variables (see Appendix B for a review.)Using (2.10), the connection between the KKU formalism and the formulation with unitary link variables can be seen rather straightforwardly, in the path-integral formulation.In the limit of infinite scalar mass m → ∞, the scalars s 1 , s 2 and s 3 are frozen to zero.Then, due to (2.10), the complex link variables x, y and z reduce to the unitary link variables up to a factor 1 . The second line of (2.1) becomes zero, and the third line becomes the magnetic plaquette term.A 'Wilson loop' on the orbifold lattice, say Tr (x ⃗ n y ⃗ n+x x⃗ n+ŷ ȳ⃗ n ), corresponds to a Wilson loop in the unitary-link formulation, say Tr U x,⃗ n U y,⃗ n+x U † x,⃗ n+ŷ U † y,⃗ n , up to an overall constant.When the scalar fields completely decouple, the flat measure for the non-compact variables reduces to the Haar measure for the compact ones.
Note that, even if the mother theory is taken to be SU(N ), the traceless condition is not imposed on each block kept by the projection.
See, however, Ref. [34] for a proposal of the removal of the U(1) part, for the Euclidean lattice.
The coupling constants for the SU(k) and U(1) parts are taken to be the same, asymptotically small value at the cutoff scale.From there, the SU(k) coupling grows toward infrared, while the U(1) coupling decreases.

Hamiltonian formalism (operator formalism)
The lattice Lagrangian obtained so far is nothing but a slightly complicated matrix model, which is just the quantum mechanics of multiple variables.The path integral is defined using the flat measure in the same way as in a harmonic oscillator, or, more generally, a system of particles in flat space.Hence, the Hamiltonian formulation can be obtained by the ordinary first-quantization procedure.We use the standard A t = 0 gauge, forcing all physical states to be gauge-invariant because of the Gausslaw constraint (2.5).
The Hamiltonian can be written in terms of the link variables x, y, z, and their canonical conjugates p x , p y , p z as follows: where (2.15) For simplicity we use xµ = x, ŷ, ẑ and pµ = px , py , pz for µ = 1, 2, 3. Then the commutation relation can be written as and (2.17) In the unitary-link formulation, the path integral is evaluated with the Haar measure, making the commutation relations more intricate.See Sec.B for details.
Another way to obtain the same result is to quantize the matrix model and then perform the orbifold projection.
Strictly speaking, some degrees of freedom cannot be gauged away.For example, when the time direction is compactified, the degrees of freedom corresponding to the Polyakov loop remain unfixed.The gauge-singlet constraint appears when these degrees of freedom are integrated out.For details, see e.g., Refs.[35,36].

Symmetry at discretized level
We now summarize the symmetries present at the discretized level, before proceeding to regularize the Hilbert space in a way suitable for representation by qubits.
• The local U(k) gauge symmetry.
Notably, the orbifold lattice construction preserves all the same symmetries as the unitary-link formulation.This alone provides us with strong motivation to study the orbifold approach further in the context of quantum simulation.
The normalization factor 1 √

Realization on a quantum computer
In this section, we discuss a realization of the pure Yang-Mills theory on the quantum computer based on the orbifold construction.
As a concrete realization, we utilize the Fock basis.Alternatively, we could consider the coordinate and momentum bases, following the treatment in Refs.[4,5].We explain the latter in Appendix C.

Fock space truncation
The Hamiltonian (2.14) does not have a free-part.We introduce two parameters µ and ω, and define 'free part' just as a trick to introduce a Fock basis: Then, we write the Hamiltonian as Note that, by definition, Ĥint is chosen as Ĥint ≡ Ĥ − Ĥfree .The 'mass' µ and 'frequency' ω are free parameters that are used to split the Hamiltonian up into free and interacting parts.In general, the efficiency of the regularization may depend on µ and ω.We define the creation and annihilation operators as and similarly for y and z.Each of x ⃗ n , y ⃗ n and z ⃗ n is a k × k complex matrix having 2k 2 real degrees of freedom, while ⃗ n labels L 3 different lattice points.Hence the number of harmonic oscillators describing the Fock space is 3 × 2k 2 × L 3 .We truncate the Hilbert space such that the excitation level of each oscillator is below Λ.The dimension of this truncated Hilbert space is Λ 6k 2 L 3 .The lattice theory is reproduced when this cutoff is removed, i.e., Λ → ∞.
Next we introduce an explicit expression in terms of qubits.We use the prescription used in Ref. [28].Let |j⟩ (j = 0, 1, • • • , Λ − 1) be the j-th excited state of the harmonic Here, we have implicitly assumed that the states under consideration are bounded both in coordinate and momentum spaces.Otherwise, this cutoff procedure may not make sense; for example, to approximate a state |x⟩ = dpe ipx |p⟩, infinitely large excitation modes are needed.This assumption is valid as long as the energy density is finite.
oscillator.We can write j in terms of binaries as j = K−1 l=0 b l 2 l .By using K ≡ log 2 Λ qubits, we can rewrite the state |j⟩ as With this encoding, the creation operator takes the form , we can express |j+1⟩⟨j| as an operator in this basis as Note that each |b ′ l ⟩⟨b l | is a linear combination of the Pauli matrices: Therefore, â † can be written as a linear combination of Pauli strings of length K = log 2 Λ (i.e. a tensor product of K Pauli spin operators).The same holds for â, and hence, x⃗ n , ŷ⃗ n and ẑ⃗ n are linear combinations of such Pauli strings.Each creation or annihilation operator consists of less than Λ 2 Pauli strings, so each four-point interaction contains at most Λ 8 Pauli strings.There are O(k 4 ) number of combinations regarding the color indices, and we must multiply this by the lattice volume (number of lattice sites) L 3 to obtain the total number of Pauli strings.Thus, the number of Pauli strings is bounded above by L 3 Λ 8 k 4 , up to a numerical constant, where each Pauli string is of length 4K = 4 log 2 Λ at most.The free part â † â can also be expressed using Pauli strings, but the cost of the free part is negligible compared to the interaction part, so we do not consider it here.

Gauge-singlet constraint
By sending the cutoff Λ to ∞, we obtain the lattice Hamiltonian acting on the extended Hilbert space containing the gauge non-singlet states.As with the Kogut-Susskind formulation, in general, it is difficult to truncate this extended Hilbert space directly to the subspace of physical gauge-invariant states.Therefore, we must choose the initial state to be a gauge singlet and simulate time-evolution precisely enough for the state to remain gauge invariant.Alternately, by adding a term like ⃗ n Tr Ĝ2 ⃗ n to the Hamiltonian, we can penalize the violation of the gauge-singlet constraint so that the gauge-singlet constraint is maintained in low-energy processes.

Ground state preparation
In the Lagrangian formulation, we imposed the condition (2.9) to obtain the desirable continuum limit.As the counterpart of this condition in the operator formulation, the gauge-invariant ground state |VAC⟩ (which is not the Fock vacuum) satisfies |VAC⟩. (3.8) Here ≃ indicates the equality up to terms that disappears in the continuum limit a → 0, as in (2.9).This does not imply, however, that x⃗ n,pq |VAC⟩ ≃ δpq |VAC⟩, as such a condition is not gauge invariant.To construct this ground state we will apply the adiabatic algorithm.The approach here will be almost identical to the construction discussed in Section 4.5 of [28], where block-encoding of the Hamiltonian is achieved from the Pauli sum form and then Wan-Kim [37] algorithm is carefully applied to efficiently prepare the ground state.For 0 ≤ s ≤ 1, we introduce Ĥ(s) as Ĥ(s) = (1 − s) Ĥfree + s Ĥ. (3.9) At s = 0 the Hamiltonian is H free defined by (3.1), and we can simply choose the Fock vacuum |0⟩ as the gauge-invariant ground state.Then, we gradually change s from 0 to 1 to prepare |VAC⟩.The quantum gate complexity of ground state preparation algorithm is analogous to analysis in Ref. [28] and given by where C ∼ L 3 k 4 Λ 8 log 2 Λ and β ∼ g 2 k 4 L 3 Λ 6 in our setup and δ is the error (in 1-norm distance) of constructed state from the true ground state.Note that this complexity is controlled by the mass gap (∆ gap ) of the adiabatic Hamiltonian H(s).Note also that, in the current setup, it is better to know the behavior of the gap in the extended Hilbert space, not just in the gauge-singlet sector, since the time evolution cannot be perfectly gauge-invariant due to various errors.One way to estimate the gap is to calculate the In this article we will not go into details, but recommend the reader to consider reading Ref. [28] For example, the truncation of the Fock space breaks the gauge invariance.Also, depending on the detail of the algorithm, the gauge invariance may be broken small amount.As mentioned before, it might be possible to avoid this issue simply by adding a term proportional to ⃗ N Tr Ĝ2 .
expectation value of the energy as a function of temperature T , by using Monte Carlo simulation of the un-gauged Euclidean theory, in which the gauge field A t is turned off.At low temperature, the energy should approach the ground state value E 0 as E(T ) = E 0 + (E 0 + ∆E)e −∆E/T + O e −(E 2 −E 0 )/T , e −2∆E/T , where ∆E is the energy gap.(Note that this calculation gives the gap at Λ = ∞.)Such analysis has already been done for the matrix model (mother theory) at s = 1 [38,39].We defer detailed consideration of non-adiabatic errors to future work.

Optimal choice of regularization parameter µ and ω
The mass µ and frequency ω in H free are parameters associated with the regularization of the Hilbert space.Depending on the choice of µ and ω, the finite-cutoff effect behaves differently as the cutoff is removed.What would be the optimal choice of µ and ω, which leads to efficient truncation?
In the orbifold construction, the continuum theory is described by the fluctuations about the background (2.9).We use the Fock states to describe this background and the fluctuations.Therefore, the wave functions of the low-lying Fock states have to be not-too-large and not-too-small, such that the background (2.9) is described efficiently.Because typical sizes of the wave function is given by ⟨n| x2 |n⟩ = n+ 1 2 µω for the n-th excited state, µω ∼ (ag 1d ) 2 is a natural choice.Then, low-lying (small-n) modes form the background efficiently, and large-n modes describe high-frequency fluctuations.

Example of efficient time-evolution algorithm
One apparent advantage of the KKU formulation is the simplicity of the Hamiltonian in the Fock basis.As explained in Sec.3.1, it takes the form Ĥ = n P.s.
where Ŝi are Pauli strings of length 4 log 2 Λ at most.The Pauli strings are unitary operators, which can easily be expressed in terms of basic quantum gates.Again, this form of the Hamiltonian is essentially the same as the one used for our matrix model paper in Ref. [28].Therefore, the same sort of algorithms can be used for efficient quantum simulations as described in Section 4.4 of Ref. [28].More specifically, one can perform block-encoding and qubitization [40] of the Hamiltonian (in Pauli sum form) and then apply Quantum Signal Processing (QSP) [41] approach.
By definition, time evolution is described by the unitary operator e −i Ĥt .The QSP uses the Jacobi-Anger expansion of the time evolution operator where J n is the Bessel function of the first kind and T n is the Chebyshev polynomial of the first kind, providing us with an efficient way to implement this expansion on a digital quantum computer.As explained in Ref. [28], the QSP, combined with the treatment of qubitization in Ref. [40], provides us with an efficient implementation of the right-hand side of (3.11).
In the implementation of the QSP, the operators R and Û constructed in Ref. [28] are multiplied to the quantum state repeatedly.In order to approximate e −i Ĥt up to error ϵ, the necessary number of applications of R and Û is O (C (||α|| • t + log ϵ −1 )) [42], where C is the cost of multiplying Û and R to a quantum state and ||α|| = n P.s.
i=1 |α i |.In our setup, C ∼ L 3 k 4 Λ 8 log 2 Λ and ||α|| ∼ g 2 1d k 4 L 3 Λ 6 in our setup.Therefore, the increase of the cost is much slower compared to the growth of the dimension of the Hilbert space, Λ 6k 2 L 3 .

Measuring glueballs
In order to perform useful computations, it is not enough to prepare a state; we must also be able to measure interesting observables.In pure Yang-Mills theory, for instance, correlation functions of glueball operators in the confining vacuum are a subject of intense interest where quantum computation may offer a distinct advantage [43][44][45].In addition, measuring the occupation numbers of glueball modes is important for computing amplitudes in gauge theory scattering processes, an exciting prospect for quantum simulation.These glueball operators are constructed as linear combinations of the Wilson loops, which are gauge-invariant, path-ordered products of unitary link variables, where the path is chosen according to the operator's representations under charge conjugation, parity, and the point group associated with the lattice.
As discussed around eq. (2.10), when the scalar fields s 1,2,3 are sufficiently suppressed by a large mass, we have If the scalar fields are not suppressed, the loop obtained from x⃗ n,µ resembles the supersymmetric Wilson loop, which is frequently considered in the context of gauge/gravity duality [46].
where U i is the unitary link variable.Thus, for a glueball operator Φ R composed of closed paths γ of length at most l, In this way, the problem of measuring correlation functions of glueball operators is reduced to computing expectation values of products of the scalar fields x ⃗ n,µ along closed loops.Note that, in this formulation, the link variables x ⃗ n,µ do not act on the Hilbert space of group elements of U(k), but are rather embedded in the larger space C k 2 , truncated appropriately to our regularization scheme.
The canonical method of measuring observables is to use the phase estimation technique introduced by Kitaev [47].This allows one to directly compute the expectation value of Φ R to precision ε using at most O(log(1/ε)) ancilla qubits and O(1/ε) controlled-U operations, where U = e iΦ R .This protocol works by applying the unitary U to the desired state |ψ⟩ a number of times depending on the state of a register of l qubits, where l = O(log(1/ε)), applying the inverse Fourier transform to the register, and conducting a measurement on the register [48].Repeating this procedure many times, one can obtain an estimate of the expectation value of the desired operator.
Alternatively, one can use the algorithm developed recently by Huang, Kueng, and Preskill [49] implementing shadow tomography to efficiently predict expectation values of local observables.This procedure carries several advantages over phase estimation; first, the algorithm uses a number of gates only polylogarithmic in the number of observables to be measured.As long as the state to be measured can be prepared efficiently, this represents an exponential speedup relative to phase estimation.Second, the protocol is well-suited for near-term quantum devices, avoiding the need for long coherence times, ancilla qubits, and controlled operations.
The algorithm predicts the expectation values of M local observables using the following method: 1. Apply a random single-qubit Clifford circuit, then measure each qubit in the computational basis.Repeat this step N = O(log(M )) times.

Use the measurement outcomes to construct N classical representations of the state using the technique of classical shadows (this can be done efficiently when
The Clifford group C n is the normalizer of the Pauli group on n qubits G n , i.e., the set of elements g which satisfy gG n = G n g.Here, a random single-qubit Clifford circuit means a tensor product of elements of C 1 over each qubit.This is equivalent to measuring each qubit in a random Pauli basis. the observables {O 1 , ..., O M } are local, see Ref. [49]), then compute the expectation value of each observable with respect to these classical representations.
3. Group the outcomes into m equal-sized groups, and compute the means of all observables for each group.For each of the M observables, the prediction for its expectation value is the median of these means.
For M observables {O 1 , ..., O M } each acting on at most r qubits, this suffices to predict {⟨O 1 ⟩, ..., ⟨O M ⟩} to precision ε with probability 1 − δ, saturating information theoretic lower bounds, where The glueball operators are composed out of Wilson loops as in Eq. (3.13).We treat the maximal length of these Wilson loops and the number of Wilson loops per glueball operator as unknown constants.In this case, the number of qubits a loop acts on is of order Kk 2 = k 2 log 2 Λ.Let f (ψ) be the time it takes to prepare the state |ψ⟩ to be measured.Then, the shadow tomography protocol suffices to measure all the specified glueball operators, {Φ R i }, to precision ε with probability 1

Other observables
Measuring glueball operators is clearly an interesting task, but we can envision numerous use cases for quantum computing where other observables are relevant.For instance, to investigate topological physics, one could introduce a theta term into the Hamiltonian and study the topological charge over the extent of a spatial lattice.For applications in nuclear physics, one may be interested in the spatial distribution of the energy or action density.These applications demand knowledge of the expectation values of a very large number of local observables, making shadow tomography especially practical.Any physical (gauge-invariant) observable can be constructed out of Wilson loops, similar to the construction of glueball operators above, so the same techniques used for measuring their expectation values may be applied.

Jordan-Lee-Preskill bound on Hilbert space
In this section, following the discussion in a series of papers by Jordan, Lee, and Preskill (JLP) [4,5], we consider how to upper-bound the probability of a state being outside the truncated Hilbert space as a function of its energy.Intuitively, this provides a justification for imposing a cutoff Λ on the local dimension of the Hilbert space, provided we do not probe physics near or above a given energy scale E(Λ).
There are two main differences between our approach and that taken by JLP.First, we truncate with respect to the Fock basis rather than the coordinate basis used by JLP.Second, in JLP, one not only has to truncate the field range, but also the number of discrete values the field is allowed to take, at least for bosonic lattice field theories.In our case, we truncate the maximal occupation number on each bosonic mode.The same arguments go through in each case with minimal modification, however.See Appendix C regarding more about the regularization in the coordinate basis.
Specifically, we truncate each oscillator at occupation number Λ, so that we require 6k 2 log 2 Λ qubits per site.We wish to show that we can simulate any physical process below an energy scale E up to error ε with Λ at most polynomial in 1/a, 1/ε, L 3 .Let P Λ be the projection operator onto the subspace where no oscillator has occupation number greater than or equal to Λ.Then, we define Here, ⃗ n labels the spatial sites, i, j label the matrix elements, ν labels the spatial direction, and σ labels the real and imaginary parts.The notation x is simply a shorthand for all of these indices.The name p out (x) refers to the probability that oscillator x is found outside the truncated Hilbert space.Let µ n(x) and σ n(x) be the mean and standard deviation, respectively, of the occuptaion number of the oscillator labeled by x.By Chebyshev's inequality, if Λ = |µ n(x) | + cσ n(x) , c > 0, then Thus, by choosing we have ⟨ψ| All we need to do, then, is upper bound max x ⟨ψ| n2 x |ψ⟩ as a function of the energy E = ⟨ψ| Ĥ |ψ⟩.This is particularly easy for Ĥfree , which is a sum of 6k 2 L 3 independent oscillators of mass µ and frequency ω.In that case, we overestimate the maximum of max x ⟨ψ| n2 x |ψ⟩ over all states |ψ⟩ with energy E by putting one oscillator in the |n⟩ state and all others in the |0⟩ state, where n = ⌈ E−ω/2 ω ⌉ (i.e., the smallest integer which ).Then, we find for the free case Λ = O kL 3/2 E √ εω .Unfortunately, in the interacting theory we find no simple analytic bounds on max x ⟨ψ| n2 x |ψ⟩ as a function of E. However, one could imagine testing this numerically on a classical or quantum computer via a well-controlled procedure.Estimating this error is crucial to understanding the resource requirements of simulating high energy scattering processes.Bounds of the type found by JLP are likely to be quite loose in practice, so there is an independent motivation for studying these errors numerically in the case of bosonic lattice field theory.
Alternatively, one could consider truncating our Hamiltonian using the coordinate and momentum basis instead of the Fock basis.In this case, it may be possible to derive upper bounds on the truncation error by bounding the fluctuations of coordinates and momenta.For an interacting theory, this bound would potentially be quite weak.In general, in order to claim stronger bounds, additional physical ingredients are needed, for instance, specifying which terms of the full Hamiltonian contribute the most energy to some physically well-motivated state.Some related comments on the efficiency of the choice of basis are given in [6].

From orbifolds to Kogut-Susskind
Both the orbifold construction (KKU) and the Kogut-Susskind formulation (KS) have pros and cons.Let us compare several aspects: • The KS requires complicated group theory associated with the harmonic expansion on the group manifold, while the KKU requires almost no group theory, using the Fock basis of ordinary harmonic oscillators instead.
This simplification in the link variables facilitates representing the Hilbert space with qubits and may open the door to the quantum simulation of lattice gauge theories with continuous variables, where harmonic oscillators are used routinely [50].
• In the KKU, the Hamiltonian is a linear combination of Pauli strings, which is easy to handle.A method of transcribing the KS formulation into Pauli operators was given by Byrnes and Yamamoto [51].However, their approach requires the calculation of Clebsch-Gordan coefficients for U(N ) and a significant amount of classical pre-processing.As a result, the form of the Pauli strings composing the U(N ) KS Hamiltonian is unclear a priori.
• No simple truncation to an orthonormal basis of the gauge-invariant Hilbert space is known in either formulation.
• The level of complication regarding the preparation of the quantum vacuum of the interacting theory appears to be the same.There may be a large difference in the efficiency of the two approaches, though, particularly with regard to the amount of truncation required for certain precision.
• In the KKU, it is straightforward to construct the U(k) theory, and it may be possible to realize SU(k), O(k), and Sp(k) theories in a similar manner.The same holds for the KS, with the reservation that the representation theory may be complicated for generic gauge groups.
• The KKU is subject to the moduli stabilization problem, while KS does not have this problem.Note, however, that in many cases, including pure Yang-Mills, the moduli stabilization problem can be resolved.
• In the KKU, some supersymmetric theories can be realized without parameter fine tuning [12].This may or may not be possible in the KS.
• The KKU may require more qubits because half of the degrees of freedom are dynamically eliminated associated with the moduli fixing.It requires careful analysis to see if this is actually the case, though.
Our conclusion is that the KKU and KS form complementary approaches.There may be problems or devices for which the orbifold approach is more natural.
At the level of the path integral, the connection between the KKU model and the unitary-link formulation is very simple, as explained in Sec.2.2.As we can see from (2.10), the complex link variables x, y and z reduce to the unitary link variables It is not trivial to determine Clebsch-Gordan coefficients.However, algebraic methods are known for SU(N ) (for small N ) and U(N ), and a numerical method is known for SU(N ) for all N [52][53][54].Quantum algorithms are also known, which appear to give speedups over classical algorithms [55,56].
(∼angular components) when the scalars (radial components) decouple.In the operator formalism, the difference between the commutation relations (2.16) and (B.5) may cause some nontrivial deviation at finite lattice spacing.
As we have seen in eq.(2.7) in Sec.2.1, the parameter g 2 1d in the orbifold lattice is related to the lattice spacing a and coupling constant g 2 4d in the 4d theory as g 2 1d = a −3 g 2 4d .Therefore, the coarse-lattice limit a → ∞ is the weak-coupling limit g 2 1d → 0. From (2.14), we can see that only the Tr|p| 2 term survives there.If we add a large scalar mass as in (2.15), the scalar part can be decoupled.Including this term, the weakcoupling limit of the orbifold lattice Hamiltonian is analogous to the strong-coupling limit of the KS formulation.
The truncation of the higher-excited modes in Fock space naturally restricts the momentum, because ⟨n| p2 |n⟩ grows linearly with n.This resembles the cutoff of the electric field (equivalently, the cutoff for the size of the representation) in the KS formulation.Both schemes naturally discard the high-energy modes in the weakcoupling limit of the KKU or the strong-coupling limit of the KS.
Finally, we mention briefly that other important reformulations of lattice gauge theory exist, notably the quantum link model approach (see the reference [10]).This kind of model might be practical when performing Monte Carlo calculations on a classical computer, where the classical lattice gauge theory action is replaced with a quantum mechanical counterpart.This is potentially useful when performing the Monte Carlo path integral in the context of quantum computation (see [57]).

Conclusion and Outlook
In this paper, we have demonstrated how to apply the orbifold methods used by KKU [12] for constructing supersymmetric lattice gauge theories to the quantum simulation of ordinary Yang-Mills theory with U(k) gauge group.In particular, we have focused on several prototypical tasks concerning the quantum simulation of high-energy processes, including preparation of the interacting vacuum state and measurement of local dynamic observables.Our construction reframes the problem of simulating the Yang-Mills theory as one of simulating a large collection of coupled harmonic oscillators.We considered the (3+1)-dimensional Yang-Mills theory for concreteness, although the generalizations to other dimensions are straightforward.We are hopeful that this novel approach will enable further exploration into the advantages of quantum computing in simulating elementary particle physics and serve as a practical formalism for some experimental quantum computing platforms.We defer a detailed comparison of the efficiency of various approaches for specific computational problems to future work and elaborate on a number of related open problems below.

Hamiltonian formulation and quantum simulation
We wish to emphasize again that the orbifold construction we have considered was originally designed for application to supersymmetric gauge theories.The Hamiltonian has similarities to supersymmetric matrix models, which are considered frequently in the context of superstring theory.Therefore, this work, together with another work by some of us [28], forms a uniform treatment of simulating high energy theories using matrix models.To simulate these models on digital quantum hardware, we can use standard Trotter methods (see an example for the Sachdev-Ye-Kitaev (SYK) model [58]) or other oracle-based algorithms (for instance, SELECT operators are widely used in quantum chemistry, qubitization, and quantum signal processing, see [59] for a comprehensive review).Alternatively, we could apply variational algorithms that are suitable for near-term quantum computers such as the variational quantum eigensolver (for an example of applications of this method to quantum field theories, see [60]).Lastly, we could perform analog Hamiltonian simulation in ultracold atomic experiments using, for instance, Rydberg atoms (see, for example, the proposal for simulating the SYK model in [61]).Thus, the Hamiltonian formulation of quantum field theories fits naturally into quantum simulation algorithms designed for quantum many-body physics and quantum chemistry.Our work provides a natural step towards solving quantum field theories using quantum devices, ie simulating high energy physics in the lab.

Towards gravity/QFT duality
Since we have discussed quantum simulation of U(k) gauge theory in this paper, it is natural to address the possibility of studying properties of large-k theory.In the context of quantum gravity, this limit is needed to approach the classical gravity regime.A natural next step is to formulate protocols for measuring correlation functions involving single, double, or multi-trace operators.Taking a proper large-k limit, one can then directly observe the approximate factorization of correlation functions, where the dual gravitational dynamics is semiclassical.Furthermore, one could study scattering problems in AdS/CFT (see, for instance, [62]) or energy spectra in gauge theories and their relation to quantum gravitational properties in the bulk (see, for instance, problems like [63]).Such studies will be especially helpful when involving supersymmetry, as in the original work by Maldacena [64].It would also be interesting to explore the QCD phase diagram in both holographic and non-holographic contexts, where way may study non-trivial phenomena such as confinement-deconfinement phase transitions, Hagedorn behaviors, and holographic Hawking-Page phase transitions [65][66][67].Direct access to the quantum states in the Hilbert space would also enable us to confirm the recently-proposed microscopic picture of the deconfinement transition [35,[68][69][70][71].It is promising that the deconfinement transition does not involve very large excitations per color degrees of freedom, as we can check via analytic calculations in weak coupling and lattice Monte Carlo simulations at strong coupling.For example, in the weak-coupling limit of the trivial vacuum of the matrix model, the average excitation level is much less than one (see e.g., Ref. [70] for explicit calculations).Therefore, a rather small value of the cutoff Λ may be enough.The topics mentioned above are important both historically and for the frontier of high energy physics, where quantum simulation could play a significant role.

Adding topological terms
Finally, we review the possibility of including topological terms.Such terms play important roles in quantum field theory, particle physics, and condensed-matter physics.In particular, so-called θ-terms exist in the standard model of particle physics and as of yet we have no satisfactory explanation for why the associated θ-angle in QCD is so small in our universe.This problem is called the strong CP problem.The leading candidate resolution is the axion scenario [72] ( see also [73]).Topological terms in Euclidean space are complex and therefore the standard approach to simulate QFTs by Markov chain Monte Carlo method suffers from the infamous sign problem.Quantum simulation, however, allows for real-time evolution, which is sign-problem-free [74].Thus, it would be interesting to extend this paper and related results to include topological terms (see a related study in classical computation [75] and a recent study about topological terms in quantum simulation [76]).The Kogut-Susskind formulation is a Hamiltonian formulation of lattice gauge theory under the A t = 0 gauge.Associated with this gauge choice, the singlet constraint is imposed on the physical states.
The Hamiltonian is given by the sum of electric and magnetic terms, Ĥ = ĤE + ĤB . (B.1) The electric term ĤE is As usual, the link variable is related to the gauge field as U µ ≃ e iagA α µ τα , where τ The plaquette can be expressed using the field strength tensor The field strength with spatial indices is the magnetic field: B 1 = F 23 , B 2 = F 31 , B 3 = F 12 .Hence, in the continuum limit a → 0, ĤE and ĤB as defined above reduce to Note that, instead of the delta function, a −3 δ ⃗ n⃗ n ′ appeared.Note also that, once (B.With this convention, it is clear that in the 'strong coupling limit' λ → ∞, the magnetic term is omitted. The operator Ûµ,⃗ n is interpreted as the coordinate of the group manifold U(k) for the link variable on the site ⃗ n in the µ-direction.Ignoring the gauge-singlet constraint, the Hilbert space is formally written as where we use the Haar measure for the integration.In other words, H µ,⃗ n = L 2 (G), where L 2 (G) is the set of square-integrable functions from G to C. Physically, this means that we consider only the normalizable states.

B.1 Realization on a quantum computer
How should we regularize this Hilbert space systematically?It would be nice if the group manifold could be discretized by a discrete subgroup, but this does not seem to work except for U(1) theory, where Z N gives an efficient discretization.
A physically elegant, but practically very hard, approach is to truncate in the loop basis.We start with the 'strong coupling limit,' where ĤB is dropped.The ground state in this limit is given by Êα µ,⃗ n |0⟩ = 0.By acting on this state with Wilson loop operators ŴC , obtained by multiplying the link variables along a closed contour C and taking their trace, for various contours, an over-complete basis of the gauge-invariant Hilbert space is obtained.The magnetic term is regarded as the smallest Wilson loop, i.e., the plaquette.When the loops do not intersect with each other or with themselves, the electric term is proportional to the sum of the lengths of the loops.When the loops intersect, the electric term joins or splits them.By identifying the Wilson loop with the string, this gives an alternative picture to the unitary link variables.A natural cutoff can be introduced by restricting the total length of a string, however, there is no known way to write down an orthonormal basis for the physical states systematically.
For quantum computation, probably the most natural option is to use the Peter-Weyl theorem, which gives the "Fourier expansion" on group manifolds: This form can be understood as follows.Above, we introduced the cutoff for the value of x, but we could use the periodic boundary condition |Λ⟩ = |0⟩ as well, assuming that the states close to the cutoff do not give non-negligible contributions.In this case, the 'shift operator' Ŝ ≡ n |n + 1⟩ ⟨n| is identified with e iδ X p, and hence, p = Ŝ1/2 − Ŝ−1/2 iδ X , up to the corrections of order δ X .From this, p2 = 2 Î− Ŝ− Ŝ−1 δ 2

X
follows.This is the same as (C.5) up to the boundary condition.Again, it is straightforward to write the right-hand side of (C.5) as a sum of Pauli strings, upon which efficient simulation algorithms may be applied.
respectively.Their bars belong to the anti-bi-fundamental representations.A ⃗ n is the gauge field for U(k) ⃗ n .See (2.3) and (2.4) for the explicit definition of the gauge transformation.
Each |j + 1⟩ ⟨j| in (3.5) is written as a sum of less than Λ = 2 K Pauli strings, because each |b ′ l ⟩ ⟨b l | in (3.6) contains one or two Pauli matrices.For example, the plaquette Tr(xyxȳ) can be written as k a,b,c,d=1 x ab y bc xcd ȳda , and hence there are k 4 combinations of (a, b, c, d).

Figure 1 .
Figure 1.Schematic picture of the embedding of the orbifold lattice in large matrices.Nonzero components (A.10) are shown in gray.