A Trailhead for Quantum Simulation of SU(3) Yang-Mills Lattice Gauge Theory in the Local Multiplet Basis

Maintaining local interactions in the quantum simulation of gauge field theories relegates most states in the Hilbert space to be unphysical -- theoretically benign, but experimentally difficult to avoid. Reformulations of the gauge fields can modify the ratio of physical to gauge-variant states often through classically preprocessing the Hilbert space and modifying the representation of the field on qubit degrees of freedom. This paper considers the implications of representing SU(3) Yang-Mills gauge theory on a lattice of irreducible representations in both a global basis of projected global quantum numbers and a local basis in which controlled-plaquette operators support efficient time evolution. Classically integrating over the internal gauge space at each vertex (e.g., color isospin and color hypercharge) significantly reduces both the qubit requirements and the dimensionality of the unphysical Hilbert space. Initiating tuning procedures that may inform future calculations at scale, the time evolution of one- and two-plaquettes are implemented on one of IBM's superconducting quantum devices, and early benchmark quantities are identified. The potential advantages of qudit environments, with either constrained 2D hexagonal or 1D nearest-neighbor internal state connectivity, are discussed for future large-scale calculations.

Precision calculations of Standard Model (SM) processes where quantum chromodynamics (QCD) [1][2][3] plays a role are important for a number of critical applications in high-energy physics (HEP) and nuclear physics (NP). These range from highly-inelastic hadron collisions and fragmentation at the energy-frontier, to neutrino-nucleus processes essential in understanding the nature of neutrinos, to searches for physics beyond the SM, to predicting the behavior of matter under extreme conditions as found in supernova and collisions of binary neutron stars. Since the discovery of QCD in the early 1970's, substantial progress has been made in developing analytic and numerical techniques to establish predictive capabilities for observables, along with a complete quantification of uncertainties, for SM processes that are beyond reach of experiment or observation. Effective field theories (EFTs) have played a central role in analytic progress, while lattice QCD [4,5] has been central to numerical progress. Generally, these techniques are frameworks for making predictions that converge to those of QCD in particular limits, and that provide systematically reducible uncertainties (both systematic and statistical), in a limited range of kinematics or SM parameters. Lattice QCD has enjoyed great success in performing precision calculations of hadronic spectra, electroweak matrix elements and form factors, properties of high-temperature low-density systems and simple multi-hadron systems (for a review of quantities, see Ref. [6]). Further, it is making inroads towards light nuclei and reactions, and higher-density systems (for recent reviews, see, for example, Refs. [7,8]). Stochastic sampling at the heart of lattice QCD calculations of finite density systems, including nuclei, encounters "sign" problems and signal-to-noise problems (also a sign problem) in computing observables. Further, real-time dynamics are not practical in such calculations (which are explicitly constructed in Euclidean space) at scale. Establishing an algorithm with efficient scaling to compute SM matrix elements with dynamical electroweak fields (chiral gauge theories) has also proven elusive. The limitations of classical computing, as demonstrated by these and other examples, were anticipated by Feynman and by Benioff in the early 1980s [9,10], and are now providing impediments to advancing important scientific applications.
Advances in the control of coherence and entanglement in the laboratory have led to the availability of first quantum devices and inspired theoretical developments leveraging the computational potential of quantum systems. Due to their natural capability of coherently manipulating quantum wavefunctions, quantum devices are expected to more naturally perform real-time dynamics of non-equilibrium and dense quantum systems at scale. One challenge in utilizing this envisioned advancement is in devising encodings of physical systems of theoretical interest onto the controllable subspace of quantum devices. As is typical, the choice of basis used to represent the structure of information impacts practical computational difficulty at subsequent steps. For the quantum simulation of lattice gauge theories, the bases that are chosen to spatially latticize the field and to digitize, or capture through finite quantum resources, the continuous local group manifold affects the processes of initialization, the complexity of the time evolution operator, creation of localized wavepackets, and the resulting distribution of information required to be captured in the final measurement process. Due to its ubiquity and potential relevance to the initialization process, the efficiency of the time evolution operator is often prioritized, and a strategy of spatially local distributions of quantum registers is employed [11,12]. With this tactic, algorithms for the implementation of large lattices can be defined through the design of a small number of local operators, commonly acting in two conjugate bases, that, when parallelized, allow for temporal propagation of the field with computation times that are independent of the volume. Thus, constructing an algorithm for the quantum simulation of fields begins with studies of the chosen basis or mapping of the field to quantum hardware.
Considerable effort has been devoted to developing quantum algorithms for the design and time evolution of lattice gauge theories on quantum devices , often through the Kogut-Susskind Hamiltonian formulation [58][59][60][61][62]. Consequently, there has been a wide range of explorations of quantum simulation basis design for fields from the scalar field to gauge theories e.g., on a position-space lattice in the eigenbasis of the field operator [11,12], in a basis of the local free-field eigenstates [63], on a lattice of momentum modes [64], in the magnetic basis [42,43], through gauge field integration in low-dimensional spaces [23,24], on an orbifold lattice [44,65], in a prepotential framework or basis of gauge invariant loop, string, and hadron excitations [32,35,48,[66][67][68][69][70][71], through the use of a spin system producing the desired continuous fields approaching a critical point [72][73][74][75][76][77][78], through discrete subgroups and group space decimation [19,25,36,79], through mesh digitization [80], using light-front formulations of lattice field theory [81,82], and in hybrid and analog approaches leveraging natural properties of trapped ions or ultracold atoms in optical lattices [15,17,18,22,33,83]. These strategies are important as optimal design is likely to depend on the physical properties of specific quantum architectures, which continue to be developed. Furthermore, this range of formulations can provide robustness in evaluating systematic uncertainties (from the performance of quantum hardware and algorithms) for observables that are inaccessible to classical computation.
In this paper, the multiplet basis utilized in the work of Byrnes and Yamamoto [13] is integrated over the local gauge space at each vertex of the lattice, reducing the Hilbert space describing the system down to the local SU(N) irreducible representations below a chosen truncation. This approach has been previously used to explore (1 + 1)-dim SU(2) lattice gauge theory [26], further implemented for a 1-dim chain of plaquettes in SU(2) lattice gauge theory [37], and is here developed for application to SU(3) lattice gauge theory.
Quantum simulations of Yang-Mills theories and QCD are in their infancy. Precision calculations of quantities that can be directly compared with experiment are far in the future, and are expected to require major advances in quantum devices, algorithms and formalism. However, in starting along the path to this ultimate objective, explorations of simple systems, establishing informative benchmarks, analyzing features of profitable mappings, observing natural structures, quantifying truncation sensitivity, and identifying amenable architectures are all important steps. We focus on understanding the behavior of simple systems, one-and two-plaquette systems, with regard to coupling, truncations in color space, the scaling of global and local basis states and operators, and the mapping of color irreps onto qubits, qutrits and qudits. We perform quantum simulations of low-truncation one-and two-plaquette systems using IBM's QExperience superconducting quantum devices. Further, we examine a framework (that appears to scale amiably) for the use of controlled-plaquette operators on qudit systems as a way to perform simulations of SU(3) Yang-Mills gauge field theory.

II. THE SU(3) YANG-MILLS HAMILTONIAN
Quantum simulations of Yang-Mills gauge theory can be performed by discretizing the gauge fields in the spatial directions using a cubic lattice of sites and defining link variables connecting adjacent sites of this underlying grid. These link variables are parallel transporters that connect, for SU (3), color vectors at one site to those at an adjacent site. The Hamiltonian is a sum over the chromo-electric and chromo-magnetic contributions, as first discussed by Kogut and Susskind [58], where g is the strong coupling constant, a is the lattice spacing between adjacent sites, and d is the number of spatial dimensions. In the irrep basis of tensor indices that are labeled by (p, q), the number of (fundamental, anti-fundamental) indices with total dimension dim(p, q) = (p + 1)(q + 1)(p + q + 2) 2 , the electric Hamiltonian is diagonal with eigenvalues determined by the Casimir operator, b |Ê (b) | 2 |p, q = p 2 + q 2 + pq + 3p + 3q 3 |p, q .
The plaquette operator,ˆ (x), is defined aŝ (x) = Tr Û 3 (x, x + aµ)Û 3 (x + aµ, x + aµ + aν)Û 3 (x + aµ + aν, x + aν)Û 3 (x + aν, x) , whereÛ 3 (x, y) are 3 × 3 unitary matrices, and µ and ν are unit vectors that define the orientation of the plaquette. In the electric basis, links are defined by states of the color irrep to which they belong, R, and the (uncorrelated) orientations in the two color spaces they connect, α and β, |R, α, β . The electric contribution from each link is proportional to the Casimir operator acting on the link without changing the color irrep, while the plaquette operators,ˆ +ˆ † , add color fluxes to the links in the plaquette, 3 and 3, which change the irrep of each link, subject to Gauss's law. Constraints imposed to define physically allowed states of the system are included through additional conditions. In the absence of external color charges and quarks, Gauss's law is satisfied by the product of link irreps at each vertex combining to a color singlet.

A. The Plaquette Operator
In the standard formulation of Hamiltonian lattice gauge theory [58], wavefunctions carry Clebsch-Gordon (CG) factors at each vertex with the effect of enforcing local gauge invariance. Using the notation of Fig. 1, an example of FIG. 1. Following an arrow convention generalizable to higher dimension, the above link labels will be employed. Indices local to one end of each link represent a set of indices characterizing the local gauge space e.g., the color spin and color hypercharge in SU (3).
the local vertex structure is (upper-left vertex) where the sum is over the quantum numbers internal to the links at the vertex. The subscript, Γ, on the SU(3) CG coefficient indexes the multiplicity of combined irreps achieved through tensor contractions. An example of this multiplicity is in the product 8 ⊗ 8 that can be combined to produce the 8-dimensional irrep in two distinct ways, symmetric and antisymmetric contractions, with distinct CGs. These multiplicities mildly complicate the calculation of plaquette matrix elements, but are otherwise benign with respect to the structure of the quantum simulation. With a truncation including only up to the single-index irreps, the vertices that contain a singlet (and are thus gauge invariant) are 1 ⊗ 1 ⊗ 1, 1 ⊗ 3 ⊗ 3, 3 ⊗ 3 ⊗ 3, and those related under global conjugation and permutation symmetries. With a truncation including the 8 irrep, described by the two index tensor with one upper and one lower index, the number of gauge invariant vertices rises to include the 1 ⊗ 8 ⊗ 8, 3 ⊗ 3 ⊗ 8 and 8 ⊗ 8 ⊗ 8.
A key role of the vertex CGs is to allow a localization of the plaquette operator, determining the magnetic Hamiltonian, as the minimal contracted loop of local link operators (directionality as in Fig. 1), U r α,β |R, a, b = ⊕R , Γ a b dim(R) dim(R ) |R , a , b R, a, r, α|R , a Γ1 R , b |R, b, r, β Γ2 , where r indicates the representation of the link operator, the a(b) label states within an irrep in the left(right) spaces, and the primes denote final state properties generated by the application of the link operator. For SU (2), the internal state labels are naturally identified with half integers capturing the third component projection of the total angular momentum. For SU (3), these labels may be the three-component vector of color isospin(T)-hypercharge(Y) rational numbers (T, T z , Y ) with T z , Y additive as utilized in Ref. [13] or, more abstractly, the Gelfand-Tsetlin patterns.
While this four-link operator is naively capable of producing transitions outside the gauge-invariant subspace, the vertex CGs prevent such transitions. To be concrete, consider the application of a plaquette operator impacting two links of a three-point vertex in an initial state of where the right shows the physical irrep configurations populated by the plaquette operator application and the associated CGs that would appear in the vertex factor. When applying theˆ † operator, 3's will be applied, according to Eq. (7), to R t and Q . Some combinations of the irreps generated by the plaquette operator are disallowed by Gauss's law, requiring information of the state of the neighboring link C 1 stored in the vertex CG to maintain gauge invariance. An example of such a configuration disallowed by the neighboring link is |3 C1 |15 Q |3 Rt as 3 ⊗ 3 does not produce a 15 or, equivalently, there is no singlet present in 3 ⊗ 3 ⊗ 15 tensor product.
As detailed in Appendix A, the vector components at each vertex can be captured analytically through the calculation of composite CG factors. As a result, the basis for quantum simulation can be simplified to expressing an SU(3) irrep on each link, leaving internal quantum numbers to impact the matrix elements comprising the local plaquette operator calculated classically. This formulation extends the observations previously made in SU(2) lattice gauge theory for a one-dimensional string of links [26] and plaquettes [37]. Defining notation through the 9-R symbol, where the sum is over all local vector and multiplicity indices, the plaquette matrix elements may be expressed as where the subscipts on the 9-R symbols graphically denote the corresponding vertex in Fig. 1. The delocalization of information necessary to consider in the application of a plaquette operator, depicted in Eq. (8), is set at the distance of neighboring links and does not grow beyond this locality for larger lattices, nor in higher dimension. Because the simple Hilbert space structure of qubit degrees of freedom will not provide the vertex CGs necessary to retain a four-link local plaquette operator, the CGs, usually separately relegated to the vertex and the operator, have been included here in their entirety. Two methods of implementation will be explored below. In Section IV, the vertex CGs of Eq. (5) will be manually captured through symmetry-projected, global wavefunctions of the lattice mapped to quantum states of a quantum device. In Section V, the vertex CGs will be captured in the structure of an eight-link operator controlled on the quantum states of the four neighboring links.

B. Connectivity in Multiplet Space
When designing operations for the implementation of dynamical processes within a Hilbert space, it is helpful to understand the natural connectivity between states. This basis-dependent feature will affect the efficiency of digital formulations of time evolution as well as their ease of implementation on quantum architectures with limited connectivity. Naturally, designing quantum hardware with connectivity matching that of the field Hilbert space (or vice versa) is expected to be advantageous.
When an SU(2) link operator in the fundamental representation acts, it is capable of raising or lowering the total angular momentum j value of the link state by ± 1 2 . When the vector components of the link Hilbert space are classically incorporated into the matrix elements of the plaquette operator, as discussed above, these j values are sufficient to describe the state of the local link degree of freedom. Thus, in a basis of multiplets, the relevant connectivity of quantum states within an SU(2) gauge link is in the form of a simple ladder, as shown in Fig. 2. While SU(2): the coefficients associated with connections between these states depend on the surrounding links and the associated local CG factors, states interact with maximally two neighboring states. For SU(3) lattice gauge theory in the multiplet basis, the connectivity among states within the local gauge link Hilbert space is only slightly more elaborate, and is well known from group theory. For the link operator in the 3 or 3, the tensor indices become, A connectivity diagram of the transitions described by Eq. (11) is shown at the bottom of Fig. 2, with black dashed lines indicating a possible truncation of irreps with up to two fundamental and two antifundamental indices. For each irrep not affected by the lower boundary of zero indices or the upper truncation, three connections exit the irrep associated with the three transitions described in Eq. (11). Connections into any irrep through the fundamental link operator also appear with maximal number three and along distinct paths from those exiting the irrep. In this sense, the one-dimensional nearest-neighbor locality of the SU(2) link-operator-generated gauge space is promoted in SU(3) to a two-dimensional hexagonal lattice in the bulk of high irrep truncation on each link. Importantly, these connections remain local upon the two-dimensional manifold.

C. Embeddings of the Gauge Space
Due to their role in defining the Hilbert space, the basis used to digitize gauge fields impacts many aspects of quantum simulations of gauge-field theories. It is for this reason that understanding the practical consequences of basis choice, or distributions of the field content onto the degrees of freedom, is expected to play a central role in optimizing simulations on different quantum architectures.
Before discussing the bases explored in this manuscript for the digitization of the SU(3) gauge field, it is worth pausing to reflect upon the basic assumption that continuous gauge fields be digitized at all. There is an alternative to digitizing the gauge field directly that retains a spatially local distribution of qubit degrees of freedom. Rather than implementing the gauge field continuum limit and then the thermodynamic and spatial continuum limits, it has been proposed, under the names of link models or qubit regularization, that this can be replaced with a one-step process by devising a spin system of appropriate local symmetries with a critical point in the same universality class as the field of interest [72][73][74][75][76][77][78]. Tuning to the lattice continuum limit at a phase transition of the latticized spin system produces an emergent relativistic field theory. One way to interpret success with this approach is through spatial blocking of the lattice producing effectively continuous field values in the continuum limit at the critical point. Just as the correlation length of the field determines the lattice volume needed to approach the thermodynamic limit for the spatial continuum with continuous fields, a second correlation length will be present in the spin lattice expressing the blocked volume necessary to capture effectively continuous fluctuations in the local field degrees of freedom. For a scalar field, it has been shown that the local digitization of the field, and the efficiency of the local quantum Fourier transform (QFT), allows the effectively continuous fluctuations in the local field degrees of freedom to be captured with double exponential convergence in local qubit number [11,12,63,84]. Benefitting from this application of the Nyquist-Shannon (NS) sampling theorem, the number of qubits per site relevant to foreseeable applications is expected to be 5 [63,84]. With current methodologies, the spatial continuum of the scalar field retains polynomial lattice artifacts as the volume-sized QFT is not expected to be efficient and thus modifications to the lattice dispersion relation will appear polynomially with the lattice spacing [63,85]. When working in a spin system, both the field continuum and the spatial continuum limits are naïvely expected to be of the latter type, allowing neither to enjoy the rapid NS convergence. It is for this reason, with an expectation that local NS convergence will retain some relevance for digitized quantum fields beyond the scalar field, that the current manuscript will focus on the local digitization of the gauge field rather than the identification of a UV completion structured as a local spin system sharing a universiality class with QCD. However, the rapidly evolving quantum ecosystem and exploratory nature of current development encourages thorough investigation of all possible avenues for embedding gauge fields into controllable quantum architectures.
Multiple embeddings of the lattice Hilbert space through the basis of SU(3) irreps will be considered in the following. The first distinction that can be made is whether the embedding is global or local. A local embedding of the Hilbert space assigns a qubit register to each link of the lattice, thus storing local information of the field in locally-distributed quantum systems. Distributing qubits/qudits across the lattice in this way is not always the most efficient use of Hilbert space. In particular there are a number of symmetries present in the gauge theory creating correlations between the link states e.g., vertex gauge symmetries, spatial parity and color-parity. Being good symmetries of the Hamiltonian, a state that begins in one symmetry sector will remain within the sector throughout its dynamical evolution. Transferring to a global embedding allows manual projections into these symmetry sectors, as utilized in Refs. [86,87]. The resulting efficient use of Hilbert space, beyond reducing the total qubit resource requirements of the calculation, protects the calculation from gauge-variant or symmetry-breaking errors. Each state of the quantum hardware is associated with a full lattice configuration with the desired projected symmetry. As long as an error maintains the computational basis, all errors maintain global and vertex symmetries of the lattice in the global embedding.
While the global embedding is advantageous for small lattices on noisy hardware, the classical computational demands for pre-conditioning the Hilbert space and computing symmetry-projected matrix elements scale poorly with the volume of the lattice. For this reason, we further explore two link-local embeddings of the lattice Hilbert space in the irrep basis. In the first approach, a single quantum register or qudit is used to capture the gauge space of each link, the operations of Section V describe how the internal modes must interact in order to express time evolution of a quantum state with respect to the magnetic Hamiltonian. While full connectivity has been permitted within a qudit, with rotations mixing populations between any two modes, it can be seen from Fig. 2 that the required mode connectivity does contain a sense of locality. In particular, the dynamical mode connectivity, reflecting that produced by the plaquette operator between irreps, will have a structure of nearest neighbor locality in a truncated and constrained two-dimensional hexagonal lattice. While this locality is an improvement over arbitrarily non-local interactions, this connectivity is potentially sub-optimal, unless a qudit architecture is designed that naturally reflects this two-dimensional structure. In particular, for a one-dimensional embedding of the irrep modes into a ladderstructured qudit, necessary mode rotations delocalize as the irrep truncation increases, reflected by the growing number of irreps per row in Fig. 2.
One natural way to address the growing two-dimensional link structure in gauge space, inspiring the second local embedding explored in this work, is to introduce a number of qudits on every link equal to the rank of the gauge group, two for SU(3), as introduced by Byrnes and Yamamoto [13]. These two qudits will reflect the tensor index structure of the irrep, in SU(3) denoted as (p, q) in Section II above, with one qudit specifying the value of p and the second indicating the value of q. In this local (p, q) representation of the local irreps, the p and q registers are simply integers from zero to the maximum number of tensor indices considered. The plaquette operator produces correlated transitions by ±1, 0 within the p and q qudits at each active link of a plaquette. In this way, further separating the link space into a pair of qudits naturally simplifies the two-dimensional connectivity shown in the left panel of Fig. 2 into a correlated pair of qudits shown in the right panel of Fig. 2, each requiring only a raising/lowering operator within a one-dimensional embedded space. With this splitting of the link space into two qudits, however, operators of up to 8 qudits controlled on another 8 qudits will be required for constructing the local time evolution operators for the one-dimensional plaquette string. As hardware-specific strategies evolve for implementing mode-isolated multi-qudit unitary rotations, tradeoffs in the fidelity of intra-and inter-qudit operations will inspire a decision: implementing a two-dimensional gauge space within a single qudit at each link versus implementing one-dimensional correlated gauge spaces within each of two qudits per link, or some spatially dependent combination of the two.
As will become clear, the choice of basis inspires different ways to perform the gauge field truncations, resulting in different convergence properties and resource requirements for simulation. Local truncations at the level of the representations in the link Hilbert spaces readily scale to larger systems, requiring resources that scale with the spatial volume of the simulation. While a global truncation connects well to intuition based on globally conserved quantities, the implementation of such a basis does not scale well with increasing system size. It would appear that an adaptive local truncation scheme constrained by a global truncation may be required for optimal use of available resources in future simulations, though constructing such a scheme is beyond the scope of this work.

III. THE SINGLE PLAQUETTE
A single plaquette is one of the simplest gauge-invariant objects that can be constructed within a lattice gauge field theory. For Yang-Mills gauge theory, the Hamiltonian responsible for its dynamics is a special case of that given in Eq. (1), given byĤ where the lattice spacing is set to a = 1, and b is an adjoint color index. The Hilbert space of a single link has been defined previously, spanned by the eigenstates of the electric-field strength operator, |R, α, β = |R, α L |R, β R . The Gauss's law constraint allows the state of the one-plaquette system to be expressed in terms of basis states of the form where R is the representation of each link. The irreducible representations of SU (3) with tensor representation T a1,a2,··· ,ap b1,b2,··· ,bq can be specified by p and q, and the global basis states for the one plaquette system is conventionally denoted by |p, q . The electric energy of a basis state is determined by the action of the Casimir operator (Eq. (3)), and is equal to four times the value of the Casimir operator on a single link multiplied by a factor of g 2 /2. The plaquette operator is defined in Eq. (4), and the subsequent section. For the one-plaquette system, matrix elements between basis states R f |ˆ |R i = 1 if R f is present in the decomposition of R i ⊗ 3, and 0 otherwise, by the completeness of CG coefficients.

A. Color Space Truncation Errors
As any computational framework is comprised of a finite number of controllable degrees of freedom, numerical explorations of gauge theories require spatial latticization as well as a form of digitization and truncation of the continuous field. As discussed above for SU(3) gauge theory on a single plaquette, digitization will here be accomplished in a gauge invariant way by truncating the number of tensor indices (p, q) at values (Λ p , Λ q ) with Λ p = Λ q a choice following the natural color parity symmetry of the system. Any truncation of the field will introduce controlled, systematic errors that must be quantified. Previous explorations of digitized scalar fields on a spatial lattice found that the field-space digitization converged double-exponentially in the number of qubits per lattice site describing the local field [11,12,63,84,88]. This convergence is attributed to the Nyquist-Shannon sampling theorem when the field and conjugate momentum bases are distributed appropriately. In Appendix B, it is shown analytically that the asymptotic form of the "color" space wavefunction for a single plaquette in SU(2) gauge theory is Gaussian with respect to the irrep dimensionality or number of tensor indices. In the case of 1 + 1 dim. SU(2) lattice gauge theory, the exponential convergence of low-lying quantities with increasing truncation has been identified previously as discussed, for example, in Ref. [48]. In this section, it is shown numerically that this Gaussian convergence is also present in the color space of an SU(3) plaquette ground state, as well as its static and dynamical observables.
Focusing upon the ground state wavefunction in a basis of irreducible representations, Fig. 3 shows an exponentially localized distribution of amplitudes. In the left panel of Fig. 3, a grid of irrep dimensionalities is established akin to that of Fig. 2. Neighboring points on this grid are connected through the magnetic plaquette operator and thus experience a sense of locality. The x-axis of this space is the difference between the number of fundamental and antifundamental indices in the tensor representation of each irrep, resulting in conjugate irreps residing at negative values and real irreps residing along the p − q = 0 axis. The dimensionalities of irreps below a truncation of Λ p = Λ q = 10 are shown explicitly, and higher index irreps would appear in the upper triangles. In this space, the ground state amplitudes of the SU(3) single plaquette wavefunction at coupling g = 0.5 are shown in the center and right panels. In the center panel, it is seen that support of this wavefunction is highly localized to the low-index regime. Of course, the extent of localization is g-dependent and becomes dispersed as g is lowered toward the weak coupling limit. The right panel presents the same wavefunction as in the center panel, but with logarithmic and quadratic functional distortions on the wavefunction amplitudes and the tensor indices, respectively. From this perspective, the asymptotic Gaussian   Fig. 2), support of the the ground state wavefunction ψ(R), shown for g = 0.5, is localized to low irrep dimensionalities (center panel). Conjugate irreps appear on the left half of the grid with real irreps appearing along the center vertical. The right panel shows log ψ(R) on a scaled quadratic grid for visual clarity of the convergence structure.    structure of the irrep-space wavefunction is visually clear. The exponential localization of the single plaquette wavefunction extends this profitable convergence also to static and dynamic observables. Figure 4 shows the convergence of the mass gap and the magnetic plaquette operator expectation value at a range of couplings. Static observables for the unit coupling are found to converge to 10 −8 percent of their asymptotic values at a low irrep truncation of Λ p = 4 up to and including tensor irreps with four fundamental and four anti-fundamental indices. As g is lowered and the wavefunction disperses in irrep space, truncation errors naturally become more dramatic. Interestingly, the mass gap demonstrates low g-dependence at high truncation, Λ p , throughout the shown coupling range. The insets of Fig. 4 provide convergence information with tensor index truncations scaled quadratically, as in the right panel of Fig. 3, such that the linear trajectories experienced at large tensor index truncations express Gaussian-type convergence structure. From these insets, one can connect necessary quantum resources to the attainable precision of local observables as the weak-coupling limit is approached. For example, percent-level precision for these quantities at couplings g ≥ 0.3 is expected to be achievable with Λ p ≤ 10 or equivalently 3-4 qubits per index register. These features are expected to apply to the link-space localization and convergence on larger lattices of SU(3) gauge theory. This suggests that SU(3) Yang-Mills simulations in a cubic spatial lattice of extent 10 × 10 × 10 could be performed with < ∼ 10 4 qubits at this coupling. It is important to keep in mind that our analysis has been performed in the electric basis, and requires increasing  resources with decreasing lattice spacing to achieve the same level of precision for any given quantity. Therefore, there is a minimum lattice spacing (coupling) below which computations are inaccessible to the electric basis for a given quantum device and available classical computing resources. Recent work by Haase et al [43] has shown in the context of QED that working instead with the magnetic basis is potentially more effective in calculations at small lattice spacings, making it an interesting area for further investigations. As a final demonstration, Fig. 5 shows similar exponential convergence properties also for dynamic observables. Calculating the time evolution of the electric energy as a function of time at increasing field truncations, Λ j , similar Gaussian precision improvements are observed. Related to the fact that the g = 1 mass gap is well captured at low truncation, the g = 1 time evolution in the left panel of Fig. 5 is well captured at low truncations even at long times. For example, the expectation value of the electric energy at time t = 100 is achievable at single precision with just 3 qubits per index register. As the coupling is reduced and the wavefunction experiences reduced locality, larger truncations are demanded to achieve precise calculations of long-time observables. The right panel of Fig. 5 quantifies this scenario for a coupling of g = 0.1. Convergence calculations such as these inform estimations of resource requirements for future lattice gauge theory simulations that will be implemented at a selection of coupling strengths and extrapolated to inform the continuum limit with a complete quantification of uncertainties. In subsequent subsections, we show the results of simulations of these systems performed on IBM's superconducting architecture.
B. Global Basis for One Plaquette

1.
8 Truncation It is informative to study a simple basis truncation of p, q ≤ 1, containing the irreps {1, 3, 3, 8}. These can be straightforwardly mapped to two qubits as {|00 , |10 , |01 , |11 } = {|1 , |3 , |3 , |8 }. In this basis, the Hamiltonian isĤ consistent with the results provided in Ref. [61] (when truncated at f 1 ≤ 2 and removing the contribution from the 6 and 6). In terms of operators acting on a two-qubit system, the electric Hamiltonian operator can be decomposed aŝ whereÎ is the identity operator. The magnetic Hamiltonian can similarly be decomposed aŝ While there is a wide range of tactics being explored for the time evolution of quantum systems [89][90][91][92][93], the method of Trotterization [94,95] is a qubit-efficient approach introducing zero auxiliary qubits. Focusing on this latter method, time evolution through Trotterization [96][97][98][99][100] for a time ∆t of a general quantum wavefunction is approximated at first and second orders as Neglecting terms proportional to the identity, the one-plaquette Hamiltonian can be separated into Trotterized operators,Ĥ =Ĥ 1 +Ĥ 2 , witĥ The matrix exponential of the first Hamiltonian contribution can be implemented with single-qubit gates, while that of the second Hamiltonian contribution can be implemented as using the decomposition of the SU(4) Cartan subalgebra [101,102].
The panels of Fig. 6 show the probability of a single plaquette remaining in the trivial vacuum, |00 , and its electric energy fluctuations for a color irrep basis truncated to {1, 3, 3, 8}. Up to four 2 nd -order Trotter steps of the form, are implemented for a coupling of g = 1. Beyond this number, it is found that the increased gate fidelity and coherence demands of the extended quantum circuit do not allow controlled mitigation of noise. The dashed curves correspond to exact classical calculations of each Trotterization, with the limit of continuous time evolution shown by the solid black curve. The data points correspond to results of the circuits discussed above implemented on IBM's Athens quantum processor [103], a superconducting qubit system in the lineage of IBM's devices using Qiskit [104]. The connectivity of this device is linear across five superconducting qubits [103] and two of the middle qubits were used to store the wavefunction of the truncated SU(3) plaquette.
The largest sources of systematic uncertainty in simulating with current quantum devices are measurements and CNOT gates, with the former dominating by a factor of ∼ 3 in basic benchmarking of the Athens device. To account for systematic errors associated with CNOT gates, previously employed extrapolation procedures [105,106] have been utilized. Mitigation of this CNOT-gate error combines the results obtained by replacing each CNOT in a circuit with an odd number, r (for r = 3, 5, 7), of CNOTs and extrapolating to r = 0 (as CNOT.CNOT=Î). Linear and quadratic extrapolations to r = 0 in the number of CNOTs-per-circuit-CNOT were performed, and 1 2 |O(linear) − O(quadratic)| was used as an estimate of the systematic uncertainty in the extrapolation of a quantity O. In many cases, the linear fit was of relatively poor quality and gate fidelity limited the reliable extraction of sufficient samples in r to estimate the systematic uncertainty at a comparison of higher polynomials. Hence, this comparison provides only an estimate and should not be considered a complete quantification of CNOT errors.
Measurement errors were mitigated in two ways. The first was implementing Qiskit's measurement filter subroutine [107] during production, which removes the leading order measurement errors by optimizing an approximate inverse of the calculated all-to-all measurement matrix. When the error introduced by application of a single CNOT gate is small compared to those of the measurement procedure, it is viable to mitigate measurement errors through the use of auxiliary qubits by implementing a majority-or unanimous-vote for the measurement result. In this democratic approach, each auxiliary qubit is connected as the CNOT target controlled on a qubit in the plaquette Hilbert space and provides one correlated measurement to inform post-selected voting. After calibration, the typical single qubit measurement error rate on the Athens processor is approximately 3% and the typical CNOT error rate is approximately 0.9% [103]. As a result, the unanimous voting criterion provides an improvement that is found in some cases to be comparable to that of the measurement filter, with degradation for circuits implemented at times distant from a calibration procedure. The initially positive results observed in this work, along with the scalability of the voting procedure, inspire future exploration of the device-dependent tuning necessary to optimize this measurement error mitigation strategy.
In addition to a choice of measurement error mitigation, the calculation shown in Fig. 6 was implemented with both a 3-and 4-CNOT gate version of e i(aX⊗X+bŶ ⊗Ŷ +cẐ⊗Ẑ) , the time evolution of the Cartan subalgebra. In the absence of noise, these two implementations should give the same results. While additional noise would be reasonably expected for the 4-CNOT calculations, temporal fluctuations in error rates of the device instead produced lower noise fluctuations for the 4-CNOT calculations. Thus, in order to express most accurately the uncertainties associated with this calculation on quantum hardware, the four implementations (3-CNOT Cartan subalgebra with unanimous voting, 4-CNOT Cartan subalgebra with the measurement filter, and others) after r-extrapolation have been combined. The uncertainty is a quadrature combination of the extrapolation errors and standard deviations of the four implementations. As a result, the uncertainties presented in Fig. 6 and throughout this manuscript are not statistical confidence intervals, but also capture the systematic errors associated with gate fidelities and temporal fluctuations of the device between calibrations that produce dominant error contributions.

6 Truncation
The next lowest-lying irreps beyond the {1, 3,3, 8} to be included in the one-plaquette basis are the 6 and6.
With six basis states, three qubits are required with two remaining unphysical states in the Hilbert space. As was leveraged in Ref. [37], the freedom of gauge-variant completion allows couplings and interactions within the unphysical subspace to be chosen to simplify the implementation of the Hamiltonian on the quantum device. The oneplaquette basis can be embedded into the Hilbert space of three qubits with the encoding, whereb = (X + iŶ )/2.
Middle qubits on the Athens quantum processor were chosen to represent the state of the system, while the two remaining qubits were used to mitigate the measurement errors of the second and fourth qubits when employing voting protocols for measurement error mitigation. A single application of the Trotterized time evolution operator acted on the trivial vacuum |000 , and the CNOT error extrapolation procedure described in Sec. III B 1 was applied. Due to the nearest neighbor couplings of the device, the Trotterized time evolution operator is decomposed into 38 CNOT gates and 37 single qubit gates. However, many of the CNOT gates in the circuit were required to compensate for the linear nearest neighbor coupling of qubits; on a device with all-to-all couplings between the qubits, this Trotterized time evolution operator could be implemented with 20 CNOT gates.
Unfortunately, implementation of this three-qubit calculation is found to exceed the capabilities of the Athens architecture with controllable systematic errors. The r ≥ 3 measurements show clear signs of coherence time saturation and the r = 1 experiences already large deviations. This combination results in an inability to perform an r extrapolation and thus an inability to mitigate the CNOT gate errors. It is possible that an extrapolation at fractional r, introducing error-exacerbating CNOT pairs at stochastically-chosen fractions of the CNOTs in the circuit as implemented in Ref. [37], could be reliably implemented, though the errors experienced already at r = 1 remain daunting anchors for extrapolation.
An additional interesting quantity to inform development is the survival probability in the physical subspace. As gauge theories are commonly designed for quantum simulation by embedding locally-interacting gauge-invariant spaces within larger Hilbert spaces, maintaining symmetry subspaces to high fidelity will be an important property of future quantum devices. These subspace fidelity demands also reside at the heart of many quantum error correction protocols with the space of logical quantum information embedded non-locally in a low-energy Hilbert space satisfying local symmetries. In Ref. [37], a gauge invariant survival probability of approximately 40% was calculated for a physical/unphysical Hilbert space ratio of 4/12 at the peak of the first oscillation in the electric energy of two plaquettes in SU(2) lattice gauge theory utilizing a circuit of 6 CNOTs on the IBM Tokyo 20-qubit quantum device. In the current application, a physical/unphysical Hilbert space ratio of 6/2 was explored with a time evolution operator of 38 CNOTs, demonstrating a survival probability of approximately 90% at the minimum of the first oscillation in the electric energy of the SU(3) global basis plaquette. While the latter represents a possible improvement in survival probability per CNOT, the former, being within both the coherence time and gate-fidelity coherence time of the device, was found to allow reliable extrapolation to a survival probability of approximately 60% and accurately capture the time evolution of electric-basis observables at the accuracy of a single Trotter step. These observations further support the necessity of multi-dimensional optimization in the design of quantum architectures. In the next subsection, the flexibility of the global basis is leveraged to perform a projection into the color parity symmetric space respected by the SU (3) Hamiltonian that is shown to allow reliable exploration of the {1, 3, 3, 8, 6, 6} dynamics through reduction onto two qubits with no unphysical subspace.

C. Color Parity Bases
To simulate the time evolution of the trivial electric vacuum, only states that are connected by repeated applications of the Hamiltonian are required to be included in the simulated basis. Observing the structure of the plaquette operator, color parity is a symmetry of SU(3) lattice Hamiltonian. Thus, time evolution will only couple states of positive color parity, to the electric strong-coupling vacuum. By including only states of the form of |R + in the basis, the evolution of the trivial electric vacuum can be simulated with a reduced Hilbert space, or a higher precision calculation can be performed using the same quantum register. 1.
The lowest non-trivial truncation in the color parity basis consists of the states |1 and |3 + = 1 √ 2 (|3 + |3 ). This can be mapped onto a single qubit with the basis choice |0 = |1 and |1 = |3 + , and the Hamiltonian becomes, With the availability of arbitrary single qubit gates, the associated time evolution can be implemented with a single unitary rotation without Trotterization. Figure 7 shows the results of performing the 3 + time evolution on the Athens quantum processor with g = 1 beginning in the electric vacuum. The combinations of measurement error and CNOT extrapolations have been employed as described in Sec. III B 1. As this calculation does not require CNOT gates, there is significantly less noise relative to the associated two-qubit calculation performed in the global basis without color parity projection. With two qubits, the color parity basis can be extended to include the |8 and |6 + states in a basis encoding of the form {|00 , |01 , |10 , |11 } = {|1 , |3 + , |6 + , |8 }, leading to a Hamiltonian of the form, To Trotterize, the single qubit terms can be grouped together, and the Cartan subalgebra (X ⊗X,Ŷ ⊗Ŷ ,Ẑ ⊗Ẑ) can be implemented as in the case of the global basis 8-truncation above. When including the 6 irrep in the color parity projected basis, there is an additionalX ⊗Ẑ term in the Hamiltonian whose Trotterized time evolution can be decomposed with the following circuit, Explicitly, the first order Trotterized time evolution operator is chosen to be implemented asÛ (t) = e −iĤ3t e −iĤ2t e −iĤ1t withĤ Implementing this Trotterized time evolution employs 10 single qubit gates and 6 CNOT gates. The accuracy of the simulation can be improved by using a second order Trotterized time evolution operator of the form, which employs 15 single qubit gates and 8 CNOT gates. Because the second order Trotter step requires fewer gates than performing two first order Trotter steps, higher order Trotterizations may be capable of improving the calculation.
Implementation of these two forms of Trotterization are presented in Fig. 8. Due to the extraX ⊗Ẑ term in the Hamiltonian, the time evolution circuit requires more CNOT gates and the additional gates in the circuit causes noise to dominate the calculation earlier than when using the global basis truncated at 8. Adding a majority choice mitigation of the measurement error in addition to the measurement filter does not significantly improve the results, indicating that the breakdown in the calculation is due to noise in the circuit rather than the measurement process. As a result, only two steps of the 1 st -order and one step of the 2 nd order Trotterizations were found to be reliable compared to the four 2 nd order steps achievable for the 8-truncated global basis. This two qubit calculation of the 6-truncated single plaquette contains all of the states coupled to the vacuum present in the three qubit global basis calculation above III B 2. However, due to the more efficient mapping, reliable time evolution is achievable with the added color parity projection.

D. Rudimentary Single Plaquette Benchmarks
While the performance of many-body dynamics cannot be captured in a single metric, benchmarks for scientific application can provide useful information toward the simulation of dynamical lattice gauge theories as quantum devices develop. Near term benchmarks are likely too rudimentary to survive into the production era, but may provide helpful guidance in the near term NISQ era. Given the exacerbated noise experienced by many quantum devices at local extrema of time evolved observables, a succinct, yet meaningful, quantity expressing device performance in this area is the extrema of the electric energy fluctuations for a single plaquette of an SU(3) lattice. Analogously to the array of hardware calibrations used to capture the high-dimensional optimization affecting the quality of operations and measurements across devices, the peaks and troughs in the fluctuation of the electric energy focuses on one informative aspect of the time evolution.
The left (right) panel of Fig. 9 shows the values of the first minimum (maximum) in the electric energy time evolution performed on the Athens quantum processor. Numerical values for the data appearing in Fig. 9 can be found in Tables V and IV of Appendix C. Being a single-qubit calculation and thus requiring no Trotterization, the data of the 3 + truncation is well controlled as seen in Fig. 7. Increasing the irreps included in the basis moving to the right also increases the number of qubits necessary to capture the Hilbert space. Within a sub-panel at fixed irrep truncation, the number of Trotter steps used to time evolve to the local extrema is increased moving to the right, increasing the gate fidelity and coherence demanded of the quantum device. Thus, from left to right each panel of Fig. 9 trades the impact of theoretical approximations for the impact of hardware noise. Ideally, this type of figure will show windows, in which Trotter errors are reduced and hardware noise has not yet overwhelmed the calculation, for an array of irrep truncations in order to inform a systematic extrapolation to the limit of infinite truncation.
As discussed in Subsection III A, the exploration of decreasing coupling increases the required Hilbert space of the ground state wavefunction in the basis of electric multiplets. While the convergence of observables at fixed g is subsequently exponential in the irrep truncation, finite computational resources, both in quantity and quality, will limit the parameter regime that can be controllably explored with extrapolation to infinite Λ irrep cut off. This  relationship has been visually translated in Fig. 9 to the presence of windows, with the smallest g reliably accessible being that for which a set of windows relevant for extrapolation are achievable. In this preliminary exploration, reliable extraction of results with increasing gauge field truncation was limited to Λ p = 1 with the {1, 3, 3, 8} basis, where Fig. 9 shows a crossing between the regime dominated by Trotter errors to the regime dominated by hardware noise with little extended intermediate window. For this reason, a coupling of g = 1 was chosen for presentation, expressing the limited Hilbert space delocalization that can be accurately captured. However, the implemented quantum circuits experience only modified rotation angles for different couplings. For this reason, it is expected that uncertainties associated with the implementation of alternate g values will be commensurate with those presented, though the angle dependence of current hardware performance supports a more thorough exploration.

E. Single Plaquette Operator Scalability
While exploring Pauli decompositions of operators are important first steps for lattice gauge theory time evolution on quantum architectures, construction of relevant operators can quickly become treacherous as the color space truncation is raised. This was experienced in the quantum simulation of SU(2) lattice gauge theory and is seen to arise also in SU (3). For this reason, exploring alternate compilation protocols for the representation of plaquette operators amenably for hardware implementation is of vital significance. In this subsection, an approach based on a decomposition into two-level unitaries will be presented for the single plaquette global basis and similar methods will be used to describe the local plaquette operator for extended lattices in Section V.
As discussed in Subsection II C, the further splitting of the local irrep basis into two registers per link representing the fundamental and antifundamental indices is likely to be practically advantageous in requiring only nearest neighbor connectivity within the two Hilbert spaces representing the link as shown at the right of Fig. 2. Because the Hilbert space of the single plaquette lattice satisfying the local Gauss's law is structurally similar to that of one (unconstrained) link of a larger lattice, the one plaquette system can be represented in a global |p, q basis with p and q digitized in a binary encoding on two separate qubit registers. With this encoding, the operators capturing the p and q index values are diagonal,p where the subscripts on Z specify which register and qubit the operator acts upon. With this representation, the electric term in the Hamiltonian becomes a sum over one-and two-qubit Pauli-Z operators as the Casimir of Eq. (3) is quadratic.
Because the connectivity of both the p-and q-registers is nearest neighbor, it is convenient to consider the (nonunitary) operator,B n that maps |p to |p − 1 and annihilates |p = |0 when p is stored in the binary encoding with n qubits. The operatorB n can be constructed recursively according tô expressed in n contributing terms with an increasing number ofb † operators in the Hilbert spaces. The operatorb has been defined previously, just below Eq. (21). Using n qubits to represent each of the p and q registers, the plaquette operator can be written asˆ where the first and second Hilbert spaces represent the p and q registers, respectively. Note again that the single plaquette lattice shares a Hilbert space structure with that of the single (unconstrained) link of an extended lattice, leading to later connections between Eq. (30) and the local link operator. As an explicit example, for a |p or |q register of four qubits, this non-unitary lowering operator would be decomposed aŝ If decomposed in the Pauli basis, theB n operator would span n k=1 2 k = 2 n+1 − 2 unique terms. Decomposing the plaquette operator in the Pauli basis subsequently demands 2 n+2 (2 n − 1) unique operators, or O(Λ p Λ q ) as the index truncation is exponential in the number of qubits per register. The Pauli decomposition of the Hermitian combination +ˆ † relevant for Trotterized time evolution presents a factor of two simplification to 2 n+1 (2 n − 1) unique operators, due to the Hermiticity of the Pauli matrices. Ignoring simplifications for the implementation of terms sharing a basis, each of these exponentially numerous terms can be implemented with O(2n) or O(log(Λ p Λ q )) CNOT gates [108] resulting in the total number of gates to implement the plaquette operator time evolution scaling exponentially with n or polynomially in Λ p,q .
While the exponential suppression of wavefunction amplitudes discussed in Subsection III A may allow this naïve approach to be practically fruitful, it is possible to restructure the plaquette time evolution decomposition for improved scaling. Again allowing p and q to be represented by two quantum registers with n qubits each, the Hermitian combination of plaquette operators present in the magnetic Hamiltonian can be written as a sum of the form whereÔ j is a tensor product operator of elements of the form {I,b,b † } ⊗2n . For example, theÔ j operators associated with the first term may consist of k identity operators, oneb operator, (n−k −1) conjugateb operators, and n identity operators for the q-register. Each of these Hermitian operators,Ô+Ô † , can be identified as a two-level unitary between computational basis states dictated by the Hilbert space locations of theb andb † operators. As discussed in Ref. [108] (4.5.2), the time evolution associated with such operators can be implemented by first transforming the basis through a Gray code [109], implementing a controlled single-qubit rotation, and then inverting the Gray code transformation. For example, implementing the time evolution associated with the last term of Eq. (31) in the |p register, contributing to the first term of Eq. (32), connects the states |1000 and |0111 in the p-register for every state in the q-register. The time evolution according to this term can then be implemented with the following circuit acting on the p-register where α is both time and coupling dependent. The second equality emphasizes the simplifications often available in the practical application of Gray code techniques, though the generic implementation of the first equality will be momentarily convenient for the scaling discussion. At the left, a Gray code is implemented in the order 0111 → 1111 → 1110 → 1100 → 1000 through the first three multi-controlled-X operators and the location of the central controlled rotation operator. With maximal Hamming distance of 2n between two bit strings spanning the {|p , |q } basis, the maximal depth of any Trotter contribution of the form e −iα(Ôj +Ô † j ) will be 2(2n) + 1 in terms of these maximally-(2n − 1)-controlledX and rotation operators. Decomposing each of these C k N OT operators into Toffoli, CNOT, and single qubit gates can be done with O(k) gates without the introduction of any auxiliary qubits [108,110]. In practice, however, the desire to avoid the demand of exponentially precise rotation gates may inspire the use of a single auxiliary qubit. This efficiency of multi-controlled CNOT operators translates directly to an equivalent efficiency in the decomposition of the general multi-controlled SU(2) rotation at the center of this circuit [111]. With this identification of two-level unitaries treated through Gray code manipulation, the total number of gates to implement the plaquette operator time evolution is found to scale polynomially with n = log 2 (Λ p + 1), the number of qubits used to represent the tensor indices of irreducible representations composing the basis.
While the qubit decompositions of these two-level contributions to the plaquette time evolution are straightforward and technically efficient, later discussions in Section V of the local plaquette operator will maintain this level of abstraction due to an expectation that currently-developing qudit frameworks may provide advantageous hardwarespecific approaches for the implementation of two-level rotations.

IV. GLOBAL BASIS: TWO PLAQUETTES
The results obtained for a single plaquette, detailed in Section III, have provided insights into the convergence of the color-representation truncation in a simple system. Some other features that are required for QCD calculations at scale only first appear in more complex systems, such as a two-plaquette system subject to spatial periodic boundary conditions (PBCs). The SU(3) two-plaquette systems are similar to those of the SU(2) system explored in Ref. [37], but with additional structure associated with the SU(3) gauge group defining the link variables. Figure 10 shows the layout of the two plaquettes, along with our conventions that define the action of plaquette operators. With an eye toward an efficient mapping of the problem onto quantum hardware, we employ the techniques used in Refs. [26,37] to "integrate over" the gauge group at each lattice site. Local gauge invariance of the theory is used to eliminate redundancies associated with the local orientations in color space, allowing the vertex amplitudes to be defined completely by the dimensionality of irreducible representations of the intersecting links. This process reduces the dimensionality of the Hilbert space and the associated resources required for quantum simulation compared with previous algorithms, for example, Ref. [13].
Similar to the methods employed for the one-plaquette system, Gauss's law can be explicitly satisfied in the global wavefunctions by construction of the basis states. Using the dimensionality of the color irrep of each link, as shown in Fig. 10, the basis states for the two-plaquette system are written as |χ(R 1 , Q 1 , R 2 , R 3 , Q 2 , R 4 ) . The gauge invariant lattice wavefunction for this two-plaquette system, as discussed in greater generality in Appendix A, is where |R, a, b is a link-state in the electric basis and R i , f, R j , k|Q k , c Γ ijk are SU(3) CG coefficients.
The global wavefunctions of the two-plaquette system are formed from combinations of these basis states, consistent with the global symmetries of the system such as: color-parity symmetry resulting from the sum of + † in the Hamiltonian, e.g., {R i , Q i } ↔ {R i , Q i }, translation invariance, and reflection symmetry. These symmetries lead to a natural block-diagonalization of the Hamiltonian in these projected bases. Quantum numbers may be assigned to the states in each block, ±1 for each of the symmetries in the case of two-plaquettes. In this section, we consider a global basis in which dynamical quantum states are mapped to symmetry-projected configurations of the full two-plaquette lattice. Two related local truncations in color space are used to explore the convergence of both local and global truncations. In limiting the local link basis to color irreps {1, 3, 3} for the two-plaquette system without constraints and symmetries, there are 3 6 independent basis states. Imposing Gauss's law at each vertex reduces this number down to 27. Further restricting to global singlet states, as is the strong coupling vacuum and preserved by the Hamiltonian, the dynamical Hilbert space becomes 9 dimensional, which decomposes into sectors of dimensions (4, 2, 2, 1) under the discrete symmetries of color parity and spatial translation. Focusing on the sector that contains the trivial vacuum, the basis states in the ++ sector are, where the superscript "++" denotes the transformation properties under color parity inversion and spatial translation, respectively.
in the −− sector. By a direct calculation of the Hamiltonian matrix elements, both the Casimir and plaquette operators, we find Hamiltonian matrices of the following form in the ++ sector, and in the other sectorsĤ The eigenvalues of these sectors are shown in the left panel of Fig. 11 as a function of the coupling, g. The axes have been re-scaled, according to their behavior in the strong and weak coupling limits, to be g 2 E vs 1/g 4 . At the left of this panel resides the strong coupling limit where the electric contributions to the Hamiltonian dominate and the ground state is well separated. At the right of this panel resides the weak coupling limit where the magnetic contributions to the Hamiltonian dominate and the ground state remains gapped below excitations. For demonstration purposes, g = 1 is chosen in what follows, however the behavior as a function of coupling should be noted when considering the lattice continuum limit, where ga → 0. While the present basis is highly truncated, and we will explore a larger basis in subsequent sections, it is orienting to see the effect of a global truncation. In the right panel of Fig. 11, the time evolution of the system initially in the trivial vacuum, |ψ (133;++) 1 for a coupling g = 1 is displayed. From this evaluation it is seen that the lowered global cutoff at a quadratic Casimir of 16 3 has an impact that increases with evolution temporal extent, a natural observation considering the low-Casimir initialization. Discrepancies first appear in magnitude at local extrema and build a significant phase shift over a few oscillations as the restricted Hilbert space of the added global truncation has effectively reduced the period. Informed by understanding of the truncation dependence of the single-plaquette wavefunction discussed in Section III, it is not surprising that the presence of states at the global truncation boundary in this system are significant.
Though severely truncated to the local {1, 3, 3} basis, this example parallels previous quantum simulations of a two plaquette SU(2) system [37] and naturally maps onto a qutrit device with each qutrit describing the color state of a link. Further, the global wavefunctions discussed above can be simulated with two qubits (embedding the four global states in the ++ sector). While current understanding indicates local bases to be advantageous at scale, global bases will continue to be valuable techniques (e.g., Refs. [86,87,112]) for exploring the quantum simulation capabilities of available quantum architectures.

Hardware Implementation
Mapping the ++ sector of Eq. (35) onto the computational basis of two qubits, |ψ A Trotterized time evolution can be constructed by further decomposing into a sum of three terms, The first order Trotterized time evolution operator used in the following implementation isÛ (t) = e −iĤ3t e −iĤ2t e −iĤ1t . Application of the first evolution contains only single qubit operators inĤ 1 , which can be implemented by single qubit quantum gates without further Trotterization, while the second evolution can be implemented using the quantum circuit in Eq. (19), and the third can be implemented with the following circuit relation The results of performing first order Trotter time steps with g = 1 beginning in the electric vacuum are shown in Fig. 12. Two middle qubits were used to store the state of the system and, when the measurement error mitigation is implemented through voting, the remaining three qubits were used to inform the post-selection described in Section III B 1. As the results show, three Trotter steps are capable of reproducing the first maximum and minimum in the evolution of the electric energy and calculations on the Athens quantum processor are in agreement with the exact calculation.   , further emphasizing a practical difference between global and local truncations. The basis states in the other sectors can be constructed straightforwardly (by inspection from the + + + states).
Using the same methods as described previously, the electric and magnetic matrix elements of the Kogut where I n is the n-dimensional identity matrix, and the full Hamiltonian in the + + + sector is the sumĤ (1338;+++) = H and higher states. Therefore, there is a fixed number of states in the global basis beyond which changes to observables from including higher-Casimir states provide an estimate of the systematic uncertainty from the local link truncation, but do not improve the fidelity of predictions. For the {1, 3, 3, 8} local truncation, the basis states 1, 2a, 2b, 3, 4, 5a, 5b have Casimirs below that of the state with the first appearance of the 6 or 6. Therefore computing observables from these wavefunctions provides a consistent prediction for the contribution from {1, 3, 3, 8} link states. Differences between predictions from the 1, 2a, 2b, 3, 4, 5a, 5b states and those of part or all of the larger Hamiltonian matrix provide an estimate of the irrep truncation uncertainties, i.e. the impact of omitting the 6, 6, 10, · · · . Figure 13 shows the time dependence of the energy in the electric field for three truncations. The dashed gray curve corresponds to the truncation imposed by the {1, 3, 3} local link truncation, the thin solid gray curve corresponds to the global truncation at a Casimir of 25/3 including the contribution from the 8 that is below the threshold for the contribution from the (6,6), and the solid black curve corresponds to the evolution from the complete matrices in Eqs. (45) and (46), locally truncated at with Λ p = Λ q = 1. The difference between the solid gray and black curves provides an estimate of the systematic uncertainty due to the truncation in color space. This parallels naive dimensional analysis that is used to estimate the systematic uncertainty introduced by the omission of counterterms in low-energy EFTs. We conclude from this analysis that the color-space truncation defined at the link level has not fully converged at g = 1, and inclusion of the 6, 6, followed by the three-index tensor representations, 10, 10, 15 and 15, will be required to obtain a result that is converged at the percent level.
A technical detail related to multiplicites in products of irreps appears in calculations with the {1, 3, 3, 8} local truncated basis, but absent in the {1, 3, 3} truncation. Specifically, as is well known, there are two distinct transitions to the 8 irrep in the tensor product 8 ⊗ 8: the symmetric and anti-symmetric contractions, with two distinct sets of CG coefficients that contribute to amplitudes, for example, in the fusion of states, A 8 +B 8 → C 8 . For the calculations in this section, the high-lying states in the spectrum involving three 8 links at one vertex, require a coherent sum over amplitudes in their fusion. Modifications to the local or global basis that would denote symmetrization or antisymmetrization at relevant vertices are not required. A consistently phased set of CG coefficients used to sum over color states at each vertex is sufficient to arrive at amplitudes and matrix elements.
In the same way that time evolution in the two-plaquette global basis truncated to {1, 3, 3} was performed above by mapping states of projected global symmetries to states in the quantum hardware, one can contemplate an analogous computation for the {1, 3, 3, 8} local truncation from the 15 × 15 matrices in Eqs. (45) and (46). The first step in developing the qubit-based quantum circuit, using the ordering of states that follows naturally from the basis of increasing global quadratic Casimir, is to project the Hamiltonian onto tensor products of Pauli and Identity operators, to give rise to coefficients of the form,Ĥ where an extra row of zeros has been added to the matrices in Eqs. (45) and (46). For the plaquette operator, there are 104 non-zero c ijkl . If the system is further truncated to eight states to be implemented on three qubits, then the number of non-zero coefficients is reduced to 30. While Pauli decompositions do not always utilize quantum resources optimally, as demonstrated in Section III E, the lack of uniform Hilbert space organization (as is present in the local (p, q) basis) leads to challenges in identifying scalable alternatives for global basis circuit decomposition. Even for this small system, involving only two plaquettes, the anticipated limitations of working with the global basis are becoming evident.

V. LOCAL BASIS: THE PLAQUETTE OPERATOR
Unlike the space-efficient global basis, where classical pre-preprocessing identifies and isolates the physical sector of the gauge field and each symmetry projected configuration of the field is mapped onto quantum hardware, the local basis distributes local qubit registers uniformly across the lattice to express local quantum numbers of the field. In this way, an operator acting on a limited number of quantum registers (dictated by its inherent spatial locality) can be developed on small lattices while subsequently retaining relevance even in the infinite volume limit. In this section, a formulation of the local plaquette operator for the SU(3) magnetic Hamiltonian is presented in the language of digital circuit elements on an architecture comprised of a qudit (d-level quantum system) representing the gauge field on each link of a one-dimensional string of plaquettes. As discussed in Section II A, the plaquette operator in a Hilbert space without naturally embedded CG factors will necessarily be controlled on the link registers neighboring the plaquette. Consider the case of a local qutrit representing the irreps {0, 1, 2} ↔ {1, 3, 3} on each link. To define the structure of the associated magnetic time evolution operator is to characterize an approach for the evolution of an infinite volume lattice with a truncation on any local excitation of the field. For each plaquette operator interacting with 8 qutrits in this system, there are 81 physical states out of the total 3 8 that satisfy Gauss's law and contain a singlet at each vertex. In order to design an operator implementing the correct quantum dynamics, it is necessary to accurately mix these physical states among themselves as well as to assure vanishing matrix elements between these states and the unphysical Hilbert space. The remaining portion of the operator, mixing the unphysical states among themselves, is a source of flexibility for the intended scientific application of gauge theory simulation and can be optimized or chosen as desired to simplify the circuit implementation. This freedom has been referred to as gauge variant completion (GVC) and will be used in the following design of the plaquette time evolution operator.
Of the 81 physical states currently being considered for plaquette operator design, there are 27 unique external link configurations, as shown in Table I. These 27 control sectors are grouped in rows by vertical and horizontal 1111  33333333  33333333  33333333  1133  1133  33113311  1313  1313  31313131  1331  1331  31313131  1333 1333 31333133 3331 331333133331 spatial parity as well as by global conjugation. The time evolution under the Hermitian plaquette operator can be implemented as a series of commuting operators in the control sectors of Table I, with a total of 27 controlled operators, mixing three physical states each, that are clearly mutually commuting.
To implement the above magnetic interaction, an architecture of qutrits is natural and amenable to generalization when higher truncations of the gauge space are designed. Consider the placement of a qutrit on each link degree of freedom. The Pauli operations in the qutrit space flip pairs of states, The natural rotation operator generalizing those available on current quantum architectures is the Givens rotation that transfers population between two levels within the qudit, Specific rotation angles of 0 and π/2 correspond to natural extensions of the Pauli-basis rotations on qubits, where the calligraphic X , Y structures are Hermitian but not unitary. In terms of these operators, one GVC of the magnetic Hamiltonian in each control sector can be constructed through the following Hermitian combinations, with coefficients and transitions determined by the plaquette matrix elements discussed is Appendix A. With the plaquette Hilbert space designated in the linearized basis of |R b |Q r |R t |Q , the active space rotations are, where, when applied in the evolution of the Yang-Mills Hamiltonian, α will be determined by g and t. As enumerated in Table I  The plaquette time evolution operator can be implemented through a collection of e −iαX X X X -type operators using the Givens rotations and a generalized CNOT operator controlling the application of a qutrit Pauli operator from Eq. (50) on the mode occupation of a second qutrit. When expressing the two-qutrit Givens rotation between the same two modes in each qutrit, a second rotation can be applied to the spectator mode to remove inadvertent rotation, where is the third mode, the complement of j, k. In order to switch the relative modes acted upon in either qudit space, a pair of X operators can be used to change basis, where X s is the X operator necessary to bring the modes of the second qutrit in line with those of the first {m, n} = {j, k} for the duration of this operator's action. This function is akin to the basis transformations commonly used in qubit implementations of multi-Pauli time evolutions. Extending the tactics of Eq. (62) to three qudits yields, where the double-circled controls represent the inclusive-or for multicontrols, applying the target operation if the state is populated in any of the controlled subspaces. For the G X X X X operator, 7 -controlled operators will be used on either side to construct the inclusive-or-controlled Pauli for removal of the rotation when the third state is populated in any of the first three qutrits where one functional realization of this multi-controlled inclusive-or operation is While this particular formulation comprised of two single-qutrit rotation operators is functionally clear and seems advantageous when considering T -costs for a potentially fault-tolerant implementation, it is entirely expected that hardware-specific variations will be made to this circuit decomposition in the course of practical implementation. Specifically, it is expected that different quantum architectures may offer unique techniques for implementing the isolated two-mode rotation when the third state is not populated in any of the four qutrits, as represents the core of this circuit decomposition. Furthermore, a qubit embedding with the common Gray-code implementation of two-level unitaries discussed in Subsection III E may be advantageous for particular architectures.
In the above, a single qudit is used to capture the local Hilbert space of each link, demanding the two-dimensional hexagonal connectivity of Fig. 2 between qudit modes. Section II C discusses the alternate opportunity to introduce two qudits per link representing the (p, q) registers of tensor indices, with the advantage that the necessary connectivity of modes in each qudit is reduced to nearest-neighbor linear, similar in form to that of a link in SU (2). The quantum circuits necessary for implementing the plaquette operator in this (p, q) basis of the local Hilbert space can be straightforwardly generalized from that above. In particular, the Hermitian combinations of Eq. (54)-(61) can be modified with the following substitutions where the {1, 3, 3} basis within a single Hilbert space is traded for the pair of Hilbert spaces |p ⊗ |q . In words, the last substitution in Eq. (67) replaces a hermitian operator mixing qudit states 0,2 (irreps 1, 3) with an identity operator in the p-register and a mixing of the lowest two levels in the q-register. The plaquette operator implemented in the (p, q) basis requires only (correlated) nearest-neighbor interactions within each of the two qudits on every link. One can write a similar translation for the control operators where the two qudit lines at the right represent the |p and |q registers from top to bottom. While the basis with a single qudit per link demonstrated homogeneity in the type of operators, G X X X X jk (t), necessary for a Trotterized implementation of the Hermitian combination of the plaquette operator and its conjugate, the (p, q) basis becomes a mix of operators G X ⊗n jk (t) for 4 ≤ n ≤ 8. These larger operators can be expressed as generalizations of Eq. (65). Explicit lists characterizing the set of operators for Trotterized plaquette implementation in both of these local bases are provided in Appendix E.

A. Comparison Between Local and Global Bases
Subsection IV A presented results for the time evolution of the two-plaquette system computed using wavefunctions defined in a global basis truncated in color space by a local link basis of {1, 3, 3}. Above, technology was developed to address this and other systems using local controlled-plaquette operators. It is valuable to assure correspondence between results of the same quantities with these two quite different approaches. Figure 14 shows the Trotterized time evolution of a |E a | 2 in the chromo-electric field starting from the trivial vacuum as a function of time with a local truncation of {1, 3, 3} for each link. The black curve corresponds to those in the right panel of Fig. 11 and Fig. 13, obtained through the time evolution of the two-plaquette system using the four global basis states in Eqs. (35) and (44). As the quantum circuits discussed above reproduce identically the unitary evolution of each component contribution to the Hamiltonian, the same curve is recovered from the local basis through Trotterized evolution leveraging the controlled-plaquette operator. For example, the last operator in Eq. (61) is implemented, schematically, as the following set of three controlled Givens rotations, and similar decompositions apply to the other contributing controlled operators. As each unitary operator is associated with a physical transition of a plaquette, with coefficients determined from matrix elements between gauge-invariant states, the Trotterization preserves gauge invariance. As usual, however, the lack of commutativity introduces discrep- for g = 1, initialized in the trivial vacuum. The black curve is calculated in the global ++ basis (the same curve as shown in the right panel of Fig. 11 and Fig. 13), while the points correspond to Trotter evolution of the contributing controlled-plaquette operators in the local basis.
ancies in the evolution, analogous to the higher dimension operators in the Symanzik action describing finite lattice spacing artifacts in lattice QCD calculations.
As anticipated, Fig. 14 shows that the results of the Trotterized evolution in the local basis converges to a welldefined function as the Trotter step size is reduced toward zero that coincides with the result from evolution using the global basis with the same imposed color truncation. This result, and others, provides a partial validation of the controlled-plaquette local basis construction for describing QCD dynamics, and thus a framework for a scalable implementation on quantum architectures.

VI. TECHNICAL ASPECTS FOR SIMULATING AT SCALE
The results of the previous sections inspire a discussion of technical hurdles that can be anticipated on the path to simulations of SU(3) lattice gauge theory at scale using local multiplet bases. We focus on issues related to scalability in volume, circuit depth, and gauge field truncation, to potential hardware implementations, and to extensions to three dimensional spatial lattices.

A. Scalability of Local Basis
This scalability discussion begins first with the total number of qubits required to express the gauge field; a diagrammatic notation parallel to that developed for large-N c scaling [113,114] is then presented to study the number of gauge invariant vertices as well as physical states and matrix elements comprising the plaquette operator. As naïvely presented above, the latter quantity can be made to directly correspond to the number of Givens rotations and thus the quantum circuit depth for the implementation of each local plaquette operator. In addition to their scaling, numerical values of these vertex and plaquette operator properties are provided to inform future practical implementations.
Qubit requirements are sensitive to the basis used to digitize and express the gauge field. The most efficient use of a hardware Hilbert space is achieved through a global basis, as discussed in Section IV, where the gauge variant space is removed through classical pre-processing and only the gauge invariant space is mapped onto quantum degrees of freedom. As demonstrated above, neither the classical pre-processing nor the subsequent compilation for quantum implementation are expected to be scalable in global bases. For this reason, local bases have been presented in Section V, trading an expanded Hilbert space for local operators that may be optimized and implemented equivalently throughout the lattice. An initial mapping of Yang-Mills to local quantum registers, as first presented by Byrnes and Yamamoto [13], requires a large number of quantum registers to define the gauge field and include the flavor, color and Dirac degrees of freedom of the quarks. Each link would be described by |p, q, T L , T z L , Y L , T R , T z R , Y R with a quantum register of qubits (or qudits) associated with each quantum number. The number of qubits for each quantum number can be determined by the cutoff, Λ p , defining a finite (p, q) space. The compression accomplished in this work, using the methods introduced in Refs. [26,37], significantly reduces the required qubit requirements for a given lattice. By integrating over the local color spaces, the number of registers per link is reduced from 8 to 2. For a lattice of L sites in each of D spatial directions, the number of qubits required for SU(3) Yang-Mills is estimated to be # qubits ∼ 2L D log 2 (Λ p + 1). For an L = 10 lattice in D = 3 dimensions truncated at Λ p = Λ q = 1 and thus restricted to color irreps {1, 3, 3, 8}, an estimate of ∼ 2000 logical qubits are required.
In the formulation of this paper, Gauss's law is implemented explicitly by the neighboring controls associated with the action of the plaquette operator between link configurations and, in particular, the dependence on these controls of the non-vanishing matrix elements that transition between plaquette configurations. Time evolution of a Yang-Mills wavefunction defined in a local basis can be accomplished by repeated applications of these controlled-plaquette operators, parallelizable at separations of two plaquettes, and thus enjoys the volume independence of the circuit depth characteristic of locally-interacting theories [96].
To estimate the scaling of quantum resources required for the implementation of each controlled-plaquette operator as the irrep truncation is raised, it is found to be convenient to work with the (p, q) coupled register mapping. With a truncation defined by the maximum number of upper and lower indices describing the highest dimension tensor within the active color space of each link, Λ p,q , it is essential to estimate the scaling as a function of increasing Λ p,q . A truncation that respects color parity, Λ p = Λ q is employed. The quasi-locality of the interactions ensures that further scaling with regard to this cutoff involves only (trivial) factors of the spacetime volume. The above analyses of the one-and two-plaquette systems indicate that, for the low-energy and low-energy-density sectors of calculations, contributions from color irreps high in the spectrum become exponentially suppressed beyond a couplingdependent value. If this suppression is maintained for extended lattices upon raising the tensor index truncation, Λ p , will dominate over the power-law scaling in the number of non-vanishing matrix elements defining the controlledplaquette operators. It is anticipated that analogous arguments will be established for localized high-energy density configurations that evolve forward in time through fragmentation and hadronization (processes that are important for nuclear and high-energy physics) to configurations of low-lying final-state hadrons. To begin to establish the scaling of the number of controlled operations that will be required for time evolution using controlled-plaquette operators acting on the local basis, it is helpful to explore and quantify the scaling of the 3-point vertex (or 2 → 1 fusion). This corresponds to multiplicities in the control of a single plaquette vertex. The maximum possible growth of the Hilbert space for this three-point vertex is Λ 6 , allowing a (Λ+1)-dimensional configuration space for each of the (p, q) registers detailing the irrep tensor structure on each link. However, not all vertex configurations satisfy Gauss's law, leading to a reduction to the number of physical vertices. Inspired by the diagrammatic tactics of large-N c scaling calculations [113,114], consider directional lines expressing propagations of fundamental or antifundamental indices. The left of Fig 15 shows the 3-point vertex characterized by p fundamental indices and q anti-fundamental indices as propagating left-right and right-left, respectively. The external lines (conventionally absent in the large-N c analysis of color-singlet objects), can be mapped directly to irreducible representations defining the local link basis. Whether or not the chosen external lines can be connected through a center region constructed by gauge invariant index contractions determines whether the vertex is physical, or contains a singlet. All gauge invariant contractions relevant to the 3-point vertex are shown at the right of Fig. 15. The SU(3) structure provides two invariant tensors: the δ-function, leading to contractions labeled c 1−6 , and the Levi-Civita, leading to contractions c 7−8 . When implementing a local truncation in the (p, q) basis, this appears as correlated constraints between the c i 's e.g., p 1 = c 1 +c 6 +c 7 ≤ Λ p . Every integer vector c refers to a physical 3-point vertex and a unique contraction pattern. However, every unique contraction pattern does not provide a unique vertex in the local (p, q) basis, as discussed in Section IV. For example, the 8 ⊗ 8 ⊗ 8 vertex can be expressed by c = {1, 0, 1, 0, 1, 0, 0, 0}, {0, 1, 0, 1, 0, 1, 0, 0}, or {0, 0, 0, 0, 0, 0, 1, 1} through a set of three δ-contractions or a pair of -contractions.
Using this diagrammatic approach, the algebraic tensorial decomposition methods of Coleman [115] indicating that the number of three-point vertex completions for given (p, q) 2,3 is or evaluations through computational packages e.g., SU-3-CG-Code Mathematica code [116,117]  irreps are given at the left of Table II. The left panel of Fig. 16 shows the results of Table II along with polynomial fits and the associated residuals. Significant reductions in fit residuals are found for polynomials up to order n = 6, beyond which the inclusion of terms of higher degree do not improve the quality of the (single precision) fit to the exact results. The results in Table II indicate Λ 6 p scaling for the number of physical 3-point vertices below a local link truncation of Λ p . This asymptotic scaling is equivalent to that of the unconstrained Hilbert space of the locally-defined vertex, indicating that satisfying Gauss's law does not modify the asymptotic polynomial scaling of the physically relevant vertices.
Another quantity whose scaling is important is the number of plaquette control sectors. This is lower bounded by the number of different products of four irreps that contains at least one singlet. The unconstrained asymptotic scaling of the associated Hilbert space is Λ 8 p . Using a similar array of techniques as those for the 3-point vertices above, explicit calculations furnish the results shown at the right of Table II    indicate that the Gauss's law constraint does, however, add an additional non-polynomial structure to the scaling, acting to reduce the physical dimensionality from the simple Λ 6,8 unconstrained values. This strong scaling provides a potentially daunting backdrop to implementation on present and future quantum devices. At low Λ p truncations, increasing the number of indices for color irreps by one in the dynamical local link basis can produce factors of ∼ 10 in the number of control sectors that are required to be implemented at each Trotter step of the plaquette time evolution.
To fully understand the circuit complexity of the local controlled plaquette operator, as discussed in Section V, with increasing gauge field truncation, the total number of non-vanishing matrix elements within the physical subspace needs to be considered. Explicit calculations of these properties are presented in Table III. While classically capturing dimensionalities for sufficiently high truncations to numerically constrain the scaling of the physical plaquette states is nontrivial, experience with the 3-and 4-point vertices suggests that an asymptotic polynomial scaling consistent with that of the unconstrained Hilbert space is likely. This would result in a scaling of Λ 16 p . Additionally, a clear lower bound of Λ 12 p can be rationalized by the freedom of the diagonal 3-point vertices before determining physically viable values for the two remaining control links. Crucially, the final step of determining the number of non-vanishing matrix elements within the physical space (connected to the time evolution circuit depth as discussed in Section V) does not contribute additional factors of Λ p to the asymptotic scaling of the number of physical plaquette states. The irrep-locality of the plaquette operator produces a surface-rather than volume-type contribution to the dimension of physical states. To see this, consider the active space of the plaquette operator contracting a 3 or 3 with the irrep at each of the four active links. As demonstrated in the connectivity diagram of Fig. 2, an application of the fundamental or anti-fundamental is capable of producing only local nearest-neighbor transitions in the structure of a two-dimensional hexagonal lattice. As such, the plaquette operator generates population in each of three new irreps for each link and thus a potential transition supported to 3 4 = 81 different final states. Many of the possible states generated will not satisfy Gauss's law at the four vertices, as enforced through the plaquette operator controls. Thus, the number of non-zero matrix elements of the controlled plaquette operator is maximally a constant factor of 81 times larger than the total number of physical states. The right column of Table III indicates that in practice, for low Λ p -truncations, this number is significantly smaller than 81, its unconstrained upper bound.
While the scaling of the number of control structures and matrix elements in the controlled-plaquette operator with cutoff in irrep space is a relatively high-order polynomial, the operator is nearly local in space, extending over just a few links. This remains the case, but involving more links, with the plethora of Hamiltonian improvements that could be implemented, for example, Refs. [85,[118][119][120]. Consequently, we expect that these operators can be determined using classical computing, and subsequently applied repeatedly throughout the lattice volume. With the anticipated color irrep localization of low-lying field configurations, we do not anticipate that classical computing resources will impose a limitation on defining the Trotterized time evolution operator for the relevant range of lattice spacing (a range that remains to be quantified). When increasing the gauge field truncation in this local multiplet basis, high-order polynomial quantum resources are traded for improvement in the physical convergence. The delocalization resulting from reducing the lattice spacing will require controlled extrapolations as the lattice scales are systematically removed. To understand these features more clearly, and to be able to better estimate resource requirements, even for modest-sized lattices, further calculations and simulations are required. For example, the use of binary encodings, as has been used for the single plaquette in Subsection III E, may lead to slower polynomial growth, a subject of future investigations.
The convergence in color space for low-lying states suggests that a low-dimension-color-irrep EFT may exist. We conjecture that plaquettes containing "high-energy links", defined by their Casimir, can be "integrated out" of the low-dimension-color-irrep space, with their effects reproduced by higher-dimension gauge-invariant operators in the Hamiltonian with coefficients determined by matching observables. This possibility remains to be explored, and will be the subject of future work. In order to perform precision calculations at scale, developing EFT techniques, such as this and those used to make predictions from lattice QCD calculations, appear to be essential.
At this point it is worth commenting on instanton configurations that mediate transitions between distinct topological sectors in Yang-Mills theory (for a review, see Ref. [121]). Far from the center of an instanton, the field strength scales as 1/r 4 from color fields that scale as 1/r 2 , and such configurations are expected to lead to modifications to the naïve quantum resource requirements. From a scaling perspective, these configurations are anticipated to introduce power-law structure in color space. The resulting convergence in the presence of an instanton in 3-dimensional calculations is expected to be simply exponential in the number of qubits, rather than a mix of single-and double-exponential convergence for the latticization and digitization, respectively, as was found in scalar field theory [11,12,63,88,122]. Furthermore, motivated by the topological freezing effects experienced when updating gauge field configurations in Euclidean lattice QCD calculations [123][124][125][126][127][128], it will be important to understand the ability of quantum simulation time evolution and state preparation techniques to efficiently capture topological charge sectors. Observed to be influential in this aspect for classical calculations, this feature further inspires the importance of thorough explorations of boundary conditions in simulation efficiency. Reliable estimates of the impact of these configurations will only become possible when 3-dim simulations become practical. However, analysis of lattice QCD gauge-field configurations, e.g., Ref. [80], in particular in regions of topological charge density may provide helpful information.
In this work, we have focused on controlled-plaquette operators for time evolution, and have not presented an explicit formulation of a state preparation. By confinement, the connected correlation functions of the vacuum are exponentially localized with a length scale set by the mass gap. As such, the techniques associated with exponentially convergent systematically-localizable operators and fixed-point quantum circuits [129,130] can, in principle, be implemented. Classical computations of a lattice system containing at least a correlation length can be used to tune the parameters of a link-based initialization quantum circuit, which can then be used throughout the larger volume of the quantum simulation. This is, of course, limited by the lattice spacing, which if small enough would exceed classical computing resources.

B. Hardware Implementation Exploratory Discussion
Focusing on the local controlled plaquette implementation due to its advantageous scaling, two mappings of the color space into the Hilbert space(s) associated with each link have been considered. The first, with a single quantum register or qudit per link, requires high connectivity among sets of 8 link registers and two-dimensional hexagonal connectivity within each. The second, with a pair of quantum registers or qudits per link, simplifies the connectivity internal to the link space from 2D-hexagonal to a correlated set of one-dimensional hierarchies requiring only (correlated) nearestneighbor raising and lowering operators within each qudit. This further organization of two registers per link, one each for the upper and lower indices of the tensor describing the color irreps of the link, technically requires additional 16register communication to implement the 8-register correlated ladder operators controlled by 8 neighboring registers.  (3) lattice gauge theory onto quad-core SRF cavities utilizing the (p, q) local basis. The light green lines indicate the lattice structure. At the left top, a one-dimensional plaquette string is illustrated in the (p, q) basis with two qudits per link. At the left bottom, the (p, q) local basis is used only for the vertical links to homogenize the quad-core operations. At the right, it is shown how an array of quad-core SRF cavities can be used to represent a two-dimensional lattice of SU(3) gauge theory in the local (p, q) basis with blue brackets indicating the cavities used to represent the (p, q) pair of qudits at local links.
While the tactics presented here can be readily implemented on qubits, qutrits, and generally qudits, they appear to be also suitable for the Superconducting Radio Frequency (SRF) cavity devices being advanced by Lawrence Livermore National Laboratory (LLNL) and Fermi National Accelerator Laboratory (FNAL) [131,132]. Relatively large cutoffs in color irrep space could be implemented for each link, even with today's cavities. It would appear that a quantum communication fabric that connects eight nearest neighbor cavities is sufficient to simulate a chain of plaquettes using controlled-plaquette operators, as detailed in this work. The SRF cavities associated with each link and those of the four external control links would be engaged coherently. This is somewhat more complicated than the standard hypercubic communication fabric used, for example, for lattice QCD calculations.
It may be profitable to develop SRF cavity systems with sufficient cavity interconnect and/or optimized control to be able to implement simulations of one-and two-plaquette SU(3) systems using the local basis (requiring 4 and 6 cavities, respectively, or twice these numbers if trading qudits for simplified intra-qudit connectivity as in the (p, q) option of the local basis), an extended one-dimensional chain of plaquettes, and a three-dimensional cube of plaquettes (requiring 12 cavities), which would provide foundational steps toward simulations of QCD. Diagrams are provided in Fig. 17 to demonstrate possible lattice connectivity through the use of quad-core SRF cavity architectures for one (left panels)-and two (right panel)-dimensional SU(3) lattices. Two different embeddings for the one-dimensional plaquette string are imagined: the first with every link represented in the (p, q) local basis with two cavities per link and the second with a mixture of the (p, q) local basis on vertical links and the single-qudit local basis on horizontal links. This mixture of local bases allows unitary operators for the magnetic time evolution to be translationally invariant among the quad-core cavities, and further emphasizes the flexibility in lattice structure amenable to hardware codesign.
Crucially, by working in localized bases, the plaquette operators designed and implemented in the process of the demonstrations suggested above remain relevant even as lattices are scaled to the infinite volume limit. Looking further into the future, one could imagine ∼ 10 3 SRF cavities with a localized communication fabric and optimized controls being used to simulate three-dimensional Yang-Mills theory on a 10 × 10 × 10 spatial lattice using or adapting techniques demonstrated in this work, and their extension to higher dimensions.
Complementarily, recent advances in forming high-fidelity qudits within the hyperfine structure of trapped-ion systems, for example Ref. [133], indicate that exploratory calculations using the qutrit encodings presented above for the {1, 3, 3} truncation may be implementable on such systems in the near future. In addition to their appealing connectivity, these systems of trapped ions are intriguingly capable of naturally generating interaction Hamiltonians that swap populations between qudit levels with a coupling constant determined by level-dependent CG sums. While the local calculation of CG coefficients is not currently foreseen to be a limiting factor, the presence of these natural interactions suggest a potential path forward for a hybrid digital-analog approach to constructing plaquette operators.
Looking forward toward future production-scale simulations, experience from classical lattice QCD calculations indicates that tuning QCD simulations, the lattice spacing and quark masses, will be required for each set of calculations. One could imagine that initial simple tunings in early productions may involve calculations of the dynamics of one-and two-plaquette systems (and higher) as benchmarks of device performance, both in execution time and fidelity.

C. Higher Spatial Dimensions
Up to this point, strings of plaquettes residing in one spatial dimension have been considered. In two spatial dimensions, four-point vertices are required, though they can be point-split into an expanded lattice of again only three-point vertices. For systems in three dimensions, vertices involve the fusion of 6 links, two in each spatial dimension. Integration over the gauge space at each vertex can be performed, and links denoted by the irrep dimensionality, as in lower dimensions. The controlled-plaquette operators will act in each of the spatial planes, controlled by the total "external" color at each vertex, as dictated by Gauss's law. Therefore, the plaquette control structure that has been developed in this work will be applicable to three dimensions, and no further structures are required. However, at each vertex, the coherent sum over all irreps that can be formed from the external four links defines the control sectors of the operator. This is reminiscent of the function provided by point-splitting, discussed extensively in the literature for such simulations, e.g., Refs. [35,134,135], here applied to the controlled plaquette operator. This is an extra layer of complexity that has to be incorporated into the quantum circuits required to analyze the system. Most of the extra layer is anticipated to be accomplished using classical computation, as it requires determining matrix elements of local objects, and does not increase in complexity with the system volume. At the level of the quantum simulation, the increase of operator structure, corresponding to an increase in the number of link projectors (and the number of orientations), the number of quantum operations scales trivially with volume.

VII. CONCLUSIONS
The quantum simulation of lattice field theories offers a path toward computing dynamical, non-equilibrium processes of importance for basic science and for advancing quantum technologies that are inaccessible to classical simulation. While classical simulations of field theories are sophisticated, with ongoing improvements to algorithms, workflows, hardware, infrastructure and community organization, quantum simulations are at their very earliest stages. Simulations of spin systems naturally map to quantum devices with quantum registers of qubits, and early real-time calculations of elastic and in-elastic processes in low-dimensions are being performed with present-day quantum devices and quantum simulators. Further, powerful formal techniques are employed, such as tensor methods, that are of benefit to both classical and quantum simulations. As three of the four fundamental forces of nature are accurately described by quantum gauge field theories-describing the interactions and dynamics of quarks, gluons, electrons, electroweak bosons, and so forth-their simulation requires the inclusion of gauge fields. Building upon a large body of work related to Hamiltonian formulations of non-Abelian quantum field theories, and somewhat recent investiga-tions related to their implementation on quantum devices, we have investigated early steps along one of the possible paths forward (a "Trailhead") for quantum simulations of lattice SU(3) Yang-Mills gauge theory, with an eye toward QCD. We anticipate that the results and insights gained in this work may be of benefit to the simulation of other lattice field theories where local symmetries play a central role and which are implemented, in part, through the action of plaquette operators. This includes algorithms that are relevant for quantum error correction for quantum hardware and computation.
In this work, the path toward implementation of QCD on quantum devices through the Kogut-Susskind Hamiltonian presented by Byrnes and Yamamoto [13] has been adapted to reduce the qubit (qudit) requirements through integration over the local gauge space at each lattice site. The continuous gauge field is digitized in the field conjugate representation by truncating the magnitude of the color-electric field that can be supported by any given link. These local discrete Hilbert spaces can be captured through a basis of quantum registers distributed locally across the lattice in two main ways: one using a mapping of irreps onto a qudit and one using two registers for the (p, q) values defining tensor indices of SU(3) irreps for each link. Applications of the plaquette operator then update each link in the plaquette, with amplitudes that depend upon irreps of the nearest neighbor links. In contrast, a global basis can be defined by mapping symmetry-projected link configurations of the entire lattice onto the states of a quantum device, and the action of the Hamiltonian determined. The resources required to implement local bases scale with the volume of the lattice, while those required for global bases scale super-polynomially. We examined, using classical simulation and IBM's superconducting quantum hardware, the dynamics and mappings of a single SU(3) plaquette onto a quantum device, and are encouraged by the exponential convergence of low-lying states and low-energy dynamics with increasing color-irrep cutoff. Two-plaquette systems were studied in detail using global and local bases and are the simplest systems that receive contributions from controlled-plaquette operators, the construction and implementation of which we have detailed. While the low-lying states and time evolution of the single-plaquette systems can be efficiently accommodated in a Hilbert space defined by qubits, the structure of the link color-irrep Hilbert space and the action of the controlled-plaquette operator lend themselves for embedding into qudit systems, such as qutrits or SRF cavity based systems.
The plaquette operator plays a central role in the simulation of lattice gauge theories, and the fidelity with which it can be implemented across the Hilbert space of a given quantum device is a key measure of the ultimate quantum simulation fidelity. Benchmarks for the performance of devices using both global bases and local bases may provide complementary information, the latter isolating the physical Hilbert space. In addition to the one-plaquette benchmarks above, analogous benchmarks for the two-, three-, four-plaquette and higher-dimensional systems will also be valuable, with the two-plaquette system providing a measure of the fidelity of the controlled-plaquette operator, the four-plaquette system sensitive to the point-split four-vertex interaction, and the three-dimensional systems sensitive to 6-link vertices. In near-term simulations, it is likely that a series of such benchmarks, starting from the single plaquette in the global basis, will be performed to identify preferred mappings onto device architectures and provide calibrations for simulations of larger systems. Furthermore, such calibrations are expected to begin the process of tuning simulation parameters, such as the cutoff in color irrep space, balancing theory approximations with device performance for the array of lattices and couplings that will eventually be necessary to extract continuum quantities.
At these early stages of development, there is value in considering, quantifying, and comparing all potential implementations of Yang-Mills gauge theories on quantum devices. Beyond the time evolution operator, the finite mass gap in Yang-Mills theories suggests that systematically localizable quantum circuits may be realizable to prepare the ground state and localized scattering states through the use of small-volume classical simulations. Unfortunately, the asymptotic polynomial scaling of the required matrix elements, while sub-exponential, is sufficiently significant to expect that large-scale classical resources will be required for both state preparation and design of the controlledplaquette operators.
Building upon previous frameworks, the multiplet basis has here been further developed through local gauge-space integrations and circuit-level decomposition of local time evolution operators, showing promise for designing scalable quantum simulations of SU(3) Yang-Mills lattice gauge field theory. Early examples of the proposed strategies have been concretely demonstrated, through implementation on a superconducting quantum architecture and through explicit enumeration of relevant operators, for one-and two-plaquette systems. These simple systems are expected to guide near-term quantum simulations, inform future codesign, and provide calibration quantities for simulations at scale.
A derivation of plaquette operator matrix elements is presented after classical incorporation of the local vector indices. Due to the single-link delocalization of the plaquette operator in a generically structured Hilbert space, as discussed in Section II A, calculating matrix elements on a three plaquette lattice will be sufficient for evaluations of plaquette strings of arbitrary length.
Below is presented the collection of physical plaquette matrix elements of Eq. (10) for the 8 unique control sectors of the local plaquette operator in the truncated qutrit space {1, 3, 3}. link in the low-lying states will have a localized distribution that is similar to Eq. (B5) up to polynomial corrections. This remains to be verified by direct simulation. Further, the comparison between SU(2) and SU (3) is complicated by the non-linear behavior in color space of SU(3). However, the expectation is that SU(3) color fields are also localized in irrep space, in a way that resembles the behavior found in SU(2).

Appendix C: Benchmarks for Single Plaquette Time Evolution
It is important to establish robust benchmarks to guide future quantum simulations of gauge field theories. Section III of the main text explores bases within the multiplet digitization of an SU(3) plaquette and their consequences for digital quantum simulation. To complement, subsection III D focuses on the effectiveness of simulations in capturing local maxima and minima in time evolution. Tables IV and V provided in this appendix show the numerical values in Fig. 9 of the main text.