SU(2) lattice gauge theory on a quantum annealer

Lattice gauge theory is an essential tool for strongly interacting non-Abelian fields, such as those in quantum chromodynamics where lattice results have been of central importance for several decades. Recent studies suggest that quantum computers could extend the reach of lattice gauge theory in dramatic ways, but the usefulness of quantum annealing hardware for lattice gauge theory has not yet been explored. In this work, we implement SU(2) pure gauge theory on a quantum annealer for lattices comprising a few plaquettes in a row with a periodic boundary condition. These plaquettes are in two spatial dimensions and calculations use the Hamiltonian formulation where time is not discretized. Numerical results are obtained from calculations on D-Wave Advantage hardware for eigenvalues, eigenvectors, vacuum expectation values, and time evolution. The success of this initial exploration indicates that the quantum annealer might become a useful hardware platform for some aspects of lattice gauge theories.


I. INTRODUCTION
Lattice gauge theory is a mainstay for studies of quantum chromodynamics (QCD) and other strongly coupled gauge theories [1]. Significant computational resources are required, but lattice calculations provide accurate and precise information about many of the interesting properties of nucleons and other hadrons directly from first principles [2]. The lattice QCD research community has a history of evaluating each type of newly available computing hardware for its possible use [3], with an emerging example being qubit-based computing which has produced a lot of enthusiasm and research activity [4]. The first simulations of a lattice gauge theory on digital qubit hardware were reported in Ref. [5] for U(1) with fermions, followed by pure gauge SU(2) [6], pure gauge SU(3) [7], and SU(2) with fermions [8]. Simulations of U(1) gauge theory have also been performed on analog quantum hardware [9][10][11]. The pure gauge SU (2) calculations are of particular relevance to the present work, along with many theoretical studies of Hamiltonian SU(2) gauge theory that relate to possible qubit approaches . No lattice gauge theory calculations on a quantum annealer have been reported until the present work, though another group has used a quantum annealer for analyzing lattice QCD results obtained from classical computers [36].
A quantum annealer is a special-purpose type of qubitbased computing device [37][38][39]. Review articles can be found in Refs. [40][41][42][43]. D-Wave Systems Inc. [44] has been building quantum annealers for several years, with each generation having a larger number of qubits and increased functionality. The current model has a quantum processing unit (QPU) with 5760 qubits, and each qubit is connected to 15 other qubits to form what is called a Pegasus architecture. Groups of physical qubits can function together as a single "logical qubit," and these logical qubits can communicate with all-to-all connectivity.
Instead of providing the user with a universal set of quantum gates, the quantum annealer is designed for a specific calculation: finding the ground state of an Ising Hamiltonian (expressed here in the quadratic unconstrained binary optimization or "QUBO" form) where each binary q i is 0 or 1 and the user can choose any real-valued coefficients h i and J ij . The hardware performs its annealing by initializing the system into the ground state of a simple Hamiltonian and then moving quasiadiabatically to the requested Ising Hamiltonian. The Ising model might seem rather far removed from the needs of lattice QCD and too restrictive for any hope of addressing a broad set of observables, but a goal of our paper is to show that this quantum annealer can indeed perform a variety of calculations for a non-Abelian lattice gauge theory. Moreover, the ability to choose directly the coefficients in Eq. (1) is a convenient alternative to what could otherwise be a long sequence of quantum gates on digital hardware. An appeal of future fault-tolerant universal quantum computers for lattice gauge theory is the potential to open avenues that appear roadblocked to classical computing methods, particularly calculating the evolution of physics in real time and calculating physics in an environment with nonzero baryon density [45]. With classical computers, standard lattice gauge theory algorithms rely on Markov chain Monte Carlo in Euclidean spacetime. The use of Euclidean spacetime offers no access to real-time dynamics. The use of Euclidean Monte Carlo means that nonzero baryon density leads to a complexvalued Monte Carlo probability distribution and therefore a sign problem [46]. Both of these roadblocks can be removed by using a Hamiltonian formulation, but then hardware requirements grow exponentially with the size of the Hilbert space. This is where the hope for a quantum computer comes in: storage of the state vectors in a qubit register can scale polynomially rather than exponentially with the system size, so combining this with an efficient quantum algorithm could lead to practical lattice gauge theory that has access to both real-time dynamics and nonzero baryon density.
Our study does not use the qubit register in a way that achieves polynomial scaling, but the flexibility to easily choose inputs to Eq. (1) does allow for significant classical preprocessing. Therefore the system size for each quantum calculation is only a portion of the physical Hilbert space, with no involvement from any unphysical Hilbert space. This is an approach that is well suited for present-day quantum annealers, leaving the goal of eliminating exponential scaling to be pursued with future hardware. Recall that a quantum annealer can be viewed as a step toward an adiabatic implementation of universal quantum computing [47][48][49] which is equivalent to the gate-based implementation [50,51], though a universal adiabatic implementation might still require more qubits than a gate-based method to achieve the equivalency. References [52,53] discuss adiabatic quantum computing for lattice gauge theory in a gate-based context.
Besides their connection to future universal quantum computers, quantum annealers should also be compared with classical computers. Classical computation will remain crucial for the development of Hamiltonian lattice gauge theory methods for years to come, and quantum annealing may be a valuable competitor for some tasks. A scaling advantage for quantum annealing relative to path integral Monte Carlo was demonstrated on D-Wave hardware in Ref. [54]. The scaling advantage for a D-Wave QPU relative to classical simulated annealing was demonstrated in Ref. [55]. For examples of speedups attainable by quantum annealing within oracular settings, see Refs. [56,57]. Error mitigation is also important for maximizing the performance of quantum annealers and this is an active research area [58]. In the present work we will show that, without any special accommodations for optimization or error mitigation, precise results can be obtained for several observables in SU(2) gauge theory on small lattices. Future studies could build on these results to learn how well quantum annealers might eventually perform relative to classical methods.
We have chosen to study the SU(2) case because it is the smallest and simplest non-Abelian gauge theory, but it is worth recalling the long-term motivations as well. SU(2) is a natural first step toward SU(3), which is the gauge group for QCD. SU (2) contributes to the understanding of SU(N ) gauge theories more generally, which helps to frame QCD in a broader perspective. SU(2) is a viable candidate for dark matter if a fermion is added [59,60].
Various aspects of SU(2) gauge theory have been studied by other researchers who are also using a Hamiltonian formulation that can connect to a qubit implementation [6,8,. Of particular interest for our work is Ref. [6], which reports the first calculation in SU(2) pure gauge theory on a quantum computer. Specifically, IBM Q Experience gate-based hardware [61] was used to compute several steps in the time evolution of an expectation value on a two-plaquette lattice. That work represents an important milestone for the community, and provides a context in which the results of our present study can be assessed.
Section II of the present work describes the formulation to be used, including the chosen truncation of gauge fields and the number of plaquettes in our lattices. The Hamiltonian matrices are constructed for these lattice systems. Section III shows how spatial symmetries can be used to block diagonalize the Hamiltonian matrices, arriving at a form that will be used as input for the quantum annealer. Section IV presents our use of the D-Wave quantum annealer for computing eigenvalues and eigenvectors as a function of the gauge coupling. Our numerical results are shown to agree with calculations from standard algorithms running on a classical computer. In Sec. V we use the D-Wave quantum annealer to calculate some vacuum expectation values as functions of the gauge coupling, compare them with classical calculations, and use them to determine the systematic effects due to gauge truncation and finite lattice size. Section VI presents a method for computing time evolution and demonstrates its performance on the D-Wave hardware. Section VII contains a summary and outlook. Appendix A contains extra information about deriving the Hamiltonian matrix, Appendix B displays some of the most important blocks from our block diagonalized Hamiltonian matrices, and Appendix C describes the adaptive quantum annealer eigensolver algorithm that we have developed for this project.

II. PREPARING THE HAMILTONIAN
The Hamiltonian for SU(2) lattice gauge theory was originally derived in Ref. [62]. Follow-up discussions in the context of quantum computing can be found in Refs. [6,14,30]. Up to an overall additive constant, the Hamiltonian iŝ whereÊ 2 i is the Casimir operator representing the chromoelectric field for the ith lattice link. We have suppressed color indices butÊ 2 i is the sum over the three (squared) color components of the standard SU(2) Lie algebra, [Ê a ,Ê b ] = i abcÊc as described, for example, in Refs. [14,30]. The plaquette operatorˆ i in Eq. (2) is the trace of the product of four gauge link operators in order (clockwise or counterclockwise) around the ith plaquette. The gauge coupling g is the only parameter in the Hamiltonian but, following Ref. [14], we have defined which will be convenient in our work. This coefficient agrees with Ref. [30] but differs from Ref. [6]. We will report energies in units of g 2 /2 so graphs of energy versus x remain bounded in the strong coupling limit, x → 0. The Hamiltonian formalism uses a spatial lattice rather than the spacetime lattices used for standard Euclidean lattice gauge theory. In the present work we employ a one-dimensional row of 2, 4, or 6 plaquettes with a periodic boundary in the long direction. Figure 1 shows the four-plaquette case.
Because SU (2) is a continuous symmetry, the Hilbert space must be truncated to encode the gauge fields into a finite qubit register. Several formulations have been considered, and Ref. [33] provides a useful comparison of the advantages and disadvantages for some of the leading options. Our work begins from the so-called angular momentum formulation that was also used in Ref. [6] though the implementation is different for a quantum annealer, as will become clear in this paper.
The color state of each gauge link can be represented by a linear combination of basis states |j, m, m where j ∈ {0, 1 2 , 1, 3 2 , . . . } identifies the irreducible representation (irrep) of SU(2) while m and m (half integers between ±j inclusive) are the SU(2) projections at each end of the link. Basis states for the entire lattice are products of these, so a basis state for Fig. 1 can be written as Color conservation (and the absence of fermions) requires that the three links arriving at any lattice site must form an SU(2) singlet, which corresponds to Gauss's law.
To apply the Hamiltonian from Eq. (2) to any basis state we need to consider both chromoelectric and plaquette terms. The chromoelectric contribution is easy to evaluate [14] and for Fig. 1 we obtain Notice that these terms are on the diagonal; a matrix element between unequal states would vanish. The plaquette contribution is a bit more involved but we provide a derivation in Appendix A. For plaquette 1 of Fig. 1, the result is where j i and J i are the SU(2) irreps in |ψ initial and |ψ final , respectively. The 6j symbols are merely square roots of ratios and are provided in Appendix A. They enforce Gauss's law automatically. Notice that |ψ initial and |ψ final will never be the same state because applying a plaquette operator necessarily changes each of those four gauge links by ± 1 2 . Therefore all plaquette terms are off diagonal.
The result in Eq. (6) applies to a lattice of any length, not just Fig. 1, because only gauge links comprising or touching the active plaquette (E,F,I,J comprise and A,B,C,D touch) are involved in the calculation. Our Eq. (6) agrees with the expression given in Ref. [6].
Because all sums over the projections m and m have already been performed to arrive at Eq. (6), any state of the lattice depends only on the irrep values j i . It is now straightforward to calculate each entry in the Hamiltonian matrix for any one-dimensional periodic lattice of plaquettes.
Step 1: Begin with the bare vacuum (all j values set to zero) and apply any number of plaquette operators to create all possible new basis states below a chosen maximum j value.
Our first example is the case considered in Ref. [6], namely the two-plaquette lattice with each gauge link truncated by j ≤ 1 2 , where the Hamiltonian matrix is The states listed to the right of the matrix identify the rows (and corresponding columns) by using the notation of Fig. 2. Our numerical studies will also go beyond this case in two ways, namely by increasing the length of the lattice and by increasing the cutoff on j. Instead of displaying these larger Hamiltonian matrices explicitly, their sizes are listed in Table I.

III. APPLYING SPATIAL SYMMETRIES
Interchanging the states 2 2 2 2 1 1 and 2 1 1 2 2 2 in Eq. (7) would leave the Hamiltonian matrix invariant. This is a The basis states of a lattice are represented in ket notation by using the multiplicities, A ≡ 2jA +1, B ≡ 2jB +1, C ≡ 2jC + 1, etc., in a pattern that matches the layout of the spatial lattice. clue to an easy block diagonalization, (8) For any matrix, performing block diagonalization is valuable because then each block can be submitted to the quantum annealer separately. This reduces the qubit requirements. It also allows the quantum annealer to provide the ground state of each block instead of only providing the ground state of the original matrix.
To generalize the block diagonalization procedure systematically to larger Hamiltonian matrices, we will apply three symmetries to the original lattice basis states: topto-bottom reflection, left-to-right reflection at a symmetry point, and spatial translation in the periodic (long) direction. It is natural to describe spatial translation by using momentum states via e ipx but recall that the coefficients in Eq. (1) need to be real, so a modified approach will be used.
As an example of applying spatial symmetries, consider the lattice with six plaquettes and j max = 1 2 . The Hamiltonian matrix is 64×64 and will not be displayed explicitly here, but we do provide the complete list of 64 basis states in Table II. These 64 states are collected into 14 sets according to spatial translation symmetry. The notation Q (j) i denotes a state in the ith set where the excitation has been translated to the right by j sites, so reading Q and jmax = 1 2 . The notation Q (j) i denotes a state in the ith set where the excitation has been translated to the right by j sites.

Set label
Number of states The starting state Notice that all 64 basis states are symmetric under a top-to-bottom reflection. That will not generally be the case for j max > 1 2 , but it is true here. If we had chosen N plaq = 2 and j max = 1 as our example then the states 1 1 3 1 1 3 and 1 3 1 1 3 1 would have been present. Block diagonalization would have been accomplished by replacing those two states with the linear combinations that are positive and negative under the reflection, i.e., The second symmetry to consider for the basis states in Table II is left-to-right reflection. Notice that most basis states are symmetric under a left-to-right reflection for some reflection point, but states in sets 7 and 8 are not. Therefore block diagonalization requires replacing pairs of those states with their symmetric and antisymmetric combinations, The third (and final) symmetry to consider is spatial translation. Invariance under spatial translation leads to conservation of linear momentum, so our Hamiltonian can be block diagonalized into blocks of definite momentum. The six allowed momenta on the six-plaquette lattice are p = −2π/3, −π/3, 0, π/3, 2π/3, and π so within the 14 sets of Table II we should multiply each state by e ipx where x is the integer location (i.e. the plaquette TABLE III. The six momenta available on a six-plaquette lattice, and the corresponding sines and cosines. number) assigned to the excitation in that state. Only differences in x really matter because any extra offset is an irrelevant overall phase. To maintain real coefficients in Eq.
(2) we can use the real and imaginary parts separately by introducing the simple factors from Table III. For the set in Table II that contains two basis states, the states to use for block diagonalizing are which correspond respectively to p = 0 and p = π. For the sets in Table II that contain three basis states, the states to use for block diagonalizing are which correspond respectively to p = 0, p = ±2π/3 real, and p = ±2π/3 imaginary. For the sets in Table II that contain six basis states, the states to use for block diagonalizing are which correspond respectively to p = 0, p = ±π/3 real, p = ±π/3 imaginary, p = ±2π/3 real, p = ±2π/3 imaginary, and p = π.
Because the p = ±π/3 and p = ±2π/3 blocks each contain forward and backward momenta, their spectra are filled with pairs of degenerate eigenvalues. We can break each of these blocks into two separate blocks by implementing a ±π/3 rotation on pairs of states, so the 18×18 block becomes a pair of 9×9 blocks, and similarly the 22×22 block becomes a pair of 11×11 blocks.
Classical computing can readily apply spatial symmetries according to the method discussed in this section. Table I lists the size of the largest block obtained after block diagonalization for the physics systems to be studied in the present work. Notice that the largest block is always the one containing the bare vacuum state 1 1 1 1 1 1 . . . and this block will be of particular interest for computing vacuum expectation values, but every block can be implemented on a quantum annealer to obtain the smallest eigenvalue and its eigenvector.

IV. COMPUTING EIGENVALUES AND EIGENVECTORS
The variational method is well known from quantum mechanics as a way to approximate the ground-state wave function and ground-state energy for a given Hamiltonian. By varying the parameters contained in the user's trial wave function, the best approximation having that particular form can be found. A more general approach is taken here. Specifically, the complete vector space (without choosing any trial wave function) will be discretized in an unbiased way and provided to the quantum annealer, which will then find the desired minimum.
In the standard variational method, the expectation value of a Hamiltonian H for any proposed state |ψ bounds the ground-state energy E 0 according to with the equality being approached as |ψ approaches the true ground state. Notice that for the overly simplistic proposal that each entry in the vector representing |ψ is either 0 or 1, we arrive at the Ising ground state of Eq. (1) which is solved by a quantum annealer. [This is true because q 2 i = q i in Eq. (1).] The general algorithm to be used in our work emerges by applying two extensions. First, a robust numerical implementation needs to handle the possibility of a null vector (q i =0 ∀ i). Second, practical applications need an implementation that can consider any proposed state |ψ without restriction to binary entries of 0 and 1. The authors of Refs. [63][64][65][66] have documented an explicit description of this general algorithm and named it the quantum annealer eigensolver (QAE). They have found the QAE to be effective for chemistry calculations on a quantum annealer, and it is also successful for lattice gauge theory as will be demonstrated presently.  (14).
The null vector is avoided by adding a penalty term as follows: where the parameter λ in the penalty term is adjusted by the user. For a small example, take H to be the 3×3 block in Eq. (8) and use three binary variables (q 1 , q 2 , q 3 ) to represent |ψ . All options for F are listed in Table IV but let us begin with the simple case of x = 0. Without the penalty term, there are two options for getting the minimum F and one of those is the unwanted null vector. Keeping the penalty term means any choice 0 < λ < 3 will provide the single normalizable state vector with the correct minimum energy. In practice there is no need to construct the explicit table because if the null vector appears then the user can scan a few λ values with the quantum annealer to find the transition point where the null vector no longer appears. The appropriate range for λ is always adjacent to that transition point. The extension beyond merely 0 and 1 entries in the state vector is accomplished by using multiple binary variables to construct a fixed-point representation. The ith entry in the proposed vector state is where K is the number of binary variables used for that entry. Notice that the a i values are evenly spaced within [−1, 1). On the quantum annealer, one logical qubit is used for each binary variable q i so finding the ground state for an N × N matrix will use N K logical qubits. Increasing K will increase the precision of the resulting eigenvalue and eigenvector.
To summarize, the state |ψ in Eq. (13) is represented by a unit vector which is (a 1 , a 2 , . . .) from Eq. (15) divided by its norm. The original Hamiltonian is used without change. See Appendix C for our implementation.
Calculations on a D-Wave quantum annealer are performed by writing python codes that call D-Wave's ocean software suite [67]. ocean provides the user with various options to explore optimizations and refinements, including the ability to adjust the annealing schedule which defines how the hardware makes the quasiadiabatic transition from its initial to final Hamiltonian. The default annealing time is 20 microseconds. The default annealing schedule during that 20 microseconds is described in Ref. [68]. We have confirmed that acceptable results are obtained for the present project by using the default time and default schedule. The only hardware parameter we need to adjust is called the chain strength.
A chain is a set of physical qubits within the D-Wave hardware that is used to represent one logical qubit. The length of each chain depends on which connections are required between this particular logical qubit and others. We allow the ocean software to automatically perform the embedding of physical qubits into logical qubits, but we must adjust the chain strength. If the strength is too low then the physical qubits within a logical qubit can disagree with one another and lead to ambiguous physics output. If the chain strength is too high then it competes with the physics terms in the intended calculation and puts a bias on the physics output. The ocean software reports every chain breaking event and this has allowed us to easily tune the chain strength to be within an acceptable range for all of our calculations.
The chain strength is implemented "behind the scenes" by the D-Wave system. Every pair of physical qubits in a chain has an implicit Hamiltonian term of the form δH = −J chain σ z j σ z k where σ z is the Pauli matrix and a subscript identifies a specific physical qubit. The coefficient J chain is the chain strength. Increasing the value of J chain increases the probability that the two qubits will be aligned [69]. For the present study, typical values are between about 1 and 5, tuned at about O(10%). Figure 3(a) shows that the ground-state eigenvalue of Eq. (8) calculated on the D-Wave Advantage quantum annealer is in agreement with classical exact diagonalization of Eq. (8). To make this graph, the coefficient of each basis state for the 3×3 block of Eq. (8) was represented by seven binary variables, so each quantum calculation used 21 logical qubits. Using somewhat fewer binary variables also gives accurate results, but the D-Wave machine has many qubits available and Fig. 3 confirms that results remain robust even when this larger number of interconnected qubits is used. Our two tunable parameters are λ and the chain strength, and approximate tuning is sufficient for each of them. The optimal range for λ is typically slightly above the eigenvalue itself so, when calculating at the sequence of x values shown in Fig. 3(a), one can use the eigenvalue obtained at one x as the initial estimate for λ at the neighboring x. The chain strength is a positive real value; for this figure we simply used one of the values 1.0, 2.0, 3.0 or 4.0 at each x location.
Each of the nine quantum annealing calculations in Fig. 3(a) used 1000 "reads" (i.e. 1000 annealing cycles) and the graph shows the smallest numerical result from each set of 1000 reads, since the smallest is always the best estimator in a variational approach. Each read used 20 microseconds of computing time on the quantum annealer. Figure 3(b) provides a histogram for the case of x = 0.5 where the chain strength was set to 2.0 and 10 of the 1000 reads had broken chains. The histogram contains the 990 unbroken cases. Because the peak of the distribution is in the bin closest to the correct energy, we would be confident of obtaining an accurate result from this quantum annealing calculation even if the classical answer had not been available.
The method used here is immediately applicable to larger physics systems (and we will do so momentarily) because each nonzero entry in any Hamiltonian matrix can be provided directly to the quantum annealer. (Zero entries never need to be provided.) This differs from the approach used in gate-based quantum computing where the Hamiltonian must be expressed as a (possibly long) sequence of quantum gates acting on a qubit register [6,33]. In addition, gate-based implementations typically include an unphysical sector in the Hilbert space that can be much larger than the physical sector [6,33] whereas our quantum annealing calculations involve only the physical Hilbert space. These quantum annealing advantages come with the notable cost of requiring many more qubits than are needed by gate-based hardware. Figure 4 displays the leading corrections arising from (a) a longer lattice and (b) a larger j max . Specifically, Fig. 4(a) shows the smallest eigenvalue from each block (determined in Sec. III) for the four-plaquette lattice with j max = 1 2 . The lowest eigenvalue comes from a 6×6 block and we use 10 4 reads per x value when running on the quantum annealer, but the other blocks are 3×3 and 1000 reads will suffice for them. The number of logical qubits per entry in the state vector is always K = 7. Figure 4(b) shows the smallest eigenvalue from each block of the two-plaquette lattice with j max = 1, again choosing K = 7. The upper eigenvalues are from a 5×5 block and two 3×3 blocks. The lowest eigenvalue comes from a 14×14 block where 10 4 reads is not enough for the original QAE algorithm [63][64][65][66] with K = 7. To continue using a maximum of 10 4 reads per calculation, we have developed an adaptive version of the QAE algorithm (we call it the AQAE algorithm) which runs first with K = 4 to find an approximate solution, then refines the solution by using K = 4 on a finer grid in the vicinity of the approximate solution, then zooms in a second time, and then a third time. This AQAE algorithm uses only K = 4 qubits per entry in the state vector but after three adaptive steps it attains the accuracy of K = 4 + 3 = 7. Additional zooms are possible until the eigenvalue is no longer improved. Data points along the lowest curve in Fig. 4(b) were obtained from the AQAE algorithm with K = 4 and using between four and nine zooms per x value. A distinction between the data points and the exact curve becomes visible on the graph for larger x values, but even the data point at x = 0.9 only deviates from the exact curve by 4%. Details about the AQAE algorithm are provided in Appendix C.
Next-order corrections are obtained by (a) extending the lattice to six plaquettes and (b) extending j max to 3 2 . Both of these cases can be studied with our same ocean code and results are displayed in Fig. 5. All of the blocks that generate the excited states in Figs. 5(a) and 5(b) are readily handled by our AQAE algorithm with 10 4 reads. The lowest block in Fig. 5(a) begins to deviate from the curve as x grows, but using more than 10 4 reads would allow that deviation to shrink. The lowest block in Fig. 5(b) is 36×36, and 10 4 reads are insufficient to see any significant improvement beyond the 14×14 results that were already shown in Fig. 4(b). Recalling that the 36×36 block contains this exact 14×14 matrix, we can conclude that the extra 22 basis states make negligible contributions at the resolution of Fig. 5(b). Therefore we show the same results from the 14×14 truncation of the ground state in Figs. 4

(b) and 5(b).
To confirm that our AQAE ocean code is working correctly, the same code was written to run on either the D-Wave quantum annealer or on a laptop with classical simulated annealing by simply changing a single integer flag in the code. We have verified that the output from classical simulated annealing is in excellent agreement with all exact curves, including the full 36×36 matrix for Fig. 5(b). Figure 6 offers an example of how the AQAE algorithm has been vital to the results obtained in this work. Without any adaptive steps, the energy eigenvalue in that graph is clearly far above the correct result when running on D-Wave hardware. There would also be a large error bar as seen by comparing the three separate calculations (each from 10 4 reads) displayed in Fig. 6. The first adaptive step provides a major improvement and successive steps continue to approach the true result, thus allowing the completion of Figs. 4 and 5 where statistical error bars are smaller than the data symbols. In contrast, the adaptive steps provide smaller improvements when calculating with a classical simulator. Figure 6 shows that classical simulated annealing from the original QAE algorithm (corresponding to no adaptive zoom steps) is already within about 5% of the correct result.

V. COMPUTING VACUUM EXPECTATION VALUES
The particles contained within SU(2) pure gauge theory are called glueballs, and their energies are obtained from differences between the eigenvalues calculated in Sec. IV, E i −E 0 , where E 0 is the smallest eigenvalue. The symmetries implemented in Sec. III identify the specific parity and momentum corresponding to each E i . Physically interesting glueball energies would be obtained from computations on larger lattices closer to the continuum limit, which is approached as the inverse gauge coupling x is increased.
The calculations in Sec. IV produced eigenvectors as well as eigenvalues, and the eigenvector corresponding to E 0 represents the theory's vacuum state. This provides access to the calculation of the vacuum-to-vacuum matrix elements that are so important in quantum field theory. In this section, vacuum expectation values are computed and used to probe the systematic effects due to lattice volume and the j max truncation.
Because we cannot use an infinite number of qubits, there is always some limit to the precision of any calculation. For the variational method, these uncertainties are O( ) for the eigenstates but O( 2 ) for the eigenvalues, where represents a perturbation. Therefore we can anticipate less precise results for vacuum expectation values than for the associated eigenvalue.
Recall from Eq. (2) that the general SU(2) Hamiltonian is the sum of a chromoelectric term and a plaquette term, where the continuum limit of the plaquette term contains a chromomagnetic contribution and an additive constant. In units of g 2 /2, the vacuum expectation value of Eq. (2) can be written as To the left of the equal sign is the smallest eigenvalue. The first(second) term on the right side can be calculated by matrix multiplication using the diagonal(off-diagonal) terms in the Hamiltonian together with the ground-state eigenvector that was computed in Sec. IV.  Figure 7 shows the three terms of Eq. (16) on a twoplaquette lattice for the available j max values, with data points obtained from the quantum annealer and curves obtained classically. Data points show small but visible deviations from the classical curves for the chromoelectric and plaquette terms separately but, as anticipated, their sum is equal to the minimum eigenvalue and is closer to its classical curve. Data points for j max = 3 2 do not appear on the graph because, as discussed in Sec. IV, those D-Wave results are not significantly resolved from the j max = 1 data points.
The full effects of gauge fields are attained as j max → ∞, and the comparison of different j max choices in Fig. 7 suggests a rapid convergence for the range of gauge couplings studied here, 0 < x < 1. The precise rate of convergence always depends on the particular observable being considered, and we see that calculations for j max = 1 and 3 2 are closer together for the full Hamiltonian than for the chromoelectric or plaquette terms separately.
To determine how results depend on lattice volume, it is convenient to divide Eq. (16) by the number of plaquettes, thereby obtaining an energy density. Classical calculations show no visible distinction between N plaq = 4 and 6 at the resolution of Fig. 8 so they appear as a single dot-dashed curve. The two-plaquette result is a nearby solid curve. Taken together, the three volumes show that these energy densities (chromoelectric, plaquette, and total) are indeed local quantities with no significant dependence on lattice volume beyond a few plaquettes, at least for the range of x considered here. As expected, computations on the D-Wave quantum annealer in Fig. 8 show smaller errors for the total energy than for the separate chromoelectric and plaquette terms.

VI. COMPUTING TIME EVOLUTION
Time evolution is a key motivation for using the Hamiltonian approach because traditional lattice gauge theory calculations employ Euclidean time and thus lack access to real-time dynamics. For quantum computing, realtime evolution can be handled with the Suzuki-Trotter approach [70,71], where e −i H is applied repeatedly for a small time step , and H denotes the Hamiltonian operator. Since a quantum annealer does not provide gates from which to build an operator, time evolution must instead be translated into a ground-state eigenvalue problem. This can be accomplished by using Kitaev-Feynman clock states [72,73] that were also used, for example, to show the equivalence of adiabatic quantum computing to gate-based quantum computing [50,51].
The basic idea is that a sequence of time values is defined, and the quantum annealer will calculate the minimum eigenvalue for the entire time sequence at once, thereby giving the state of the system at all times. A clear derivation can be found in Ref. [74] where the approach is called the time-embedded discrete variational principle (TEDVP). Reference [74] shows that the functional to be minimized is where |Ψ t is the state of the system at time t, |t is the state of the clock at time t, and the clock Hamiltonian is where U t = e −i Ht performs the evolution from t to t + and C 0 is a penalty term used to specify the initial state at time t = 0. Because Eq. (17) has the same form as Eq. (14) up to an additive constant, we can apply the QAE directly to the TEDVP calculation of time evolution. Notice that the clock Hamiltonian acts on a compound state built from both the physical state and the clock state. The compound state is larger than the physical state by a factor of the number of time steps, so the number of qubits required will increase by this same multiplicative factor. Implementing a QAE+TEDVP algorithm on the quantum annealer is hampered by imaginary terms in U t because the D-Wave hardware needs real entries throughout Eq. (1). We will handle this by working in a basis where the Hamiltonian is purely imaginary.
Consider the case of a two-plaquette lattice that is described sufficiently accurately by the truncated Hamiltonian shown in Eq. (8). The 3×3 block describes the ground state plus two excited states, and we want to calculate the oscillation between those two excited states as a function of time. As indicated by the labels in Eq. (8), this will be oscillation between a superposition of singleplaquette excitations 1 √ 2 2 2 2 2 1 1 + 2 1 1 2 2 2 and a pair of round-the-world excitations 1 2 2 1 2 2 . Both options have exactly four excited gauge links, so they are exactly degenerate at x = 0 and nearly degenerate for small x.
There is no change of basis that converts this 3×3 block into a 3×3 imaginary matrix. However, in the strong coupling (small x) region where Eq. (8) applies, we can augment it with an additional heavy state to form this matrix, that can be written (up to a constant) in a purely imaginary form, h ± = 1 2 18 + 33x 2 ± 65x 4 + 1116x 2 + 324 . Since our calculation of oscillations will only depend on energy differences, the constant can be subtracted from Eq. (20), leaving an imaginary matrix that is block diagonal. The block with h + contains the ground state and the fictitious heavy state. The block with h − contains the two intermediate states that are of interest to us, and its time evolution is represented by where ω = h − g 2 /2. Note that, for small x, the effect of adding the fictitious heavy state [which has the value 6 on the diagonal of Eq. (19)] is smaller than the effect of truncating the Hamiltonian matrix down to its 3×3 form [because Eq. (B3) has values smaller than 6 on the diagonal that must be truncated to arrive at the 3×3 form]. Therefore, since we are considering an example where the 3×3 truncation is sufficiently accurate, we are justified to add the heavy state.
It is now straightforward to implement our QAE+TEDVP code on D-Wave hardware.
However, today's quantum annealers are only intended to handle situations that do not have sign problems. Specifically, D-Wave hardware is designed to handle "stoquastic" Hamiltonian matrices, which have only nonpositive off-diagonal elements in the computational basis [69], but Eq. (18) does not have this form when U t is given by Eq. (22) with an arbitrary time step . Figure 9 presents the output from several runs of our code that used only two times each: t = 0 and a larger t = .
With two states at two times and K = 7, each computation used 28 logical qubits. To obtain the precision of Fig. 9, we used 5 × 10 4 reads. However, even noisy data would be sufficient to provide a useful estimate of the frequency ω directly from the D-Wave data. This is valuable because we can then choose the truly stoquastic case of = π/ω for followup computations, where D-Wave hardware will handle time evolution well. Figure 10 shows 12 time values obtained from a single run of stoquastic time evolution using our QAE+TEDVP algorithm on the D-Wave Advantage hardware. This computation used K = 2 for each of the two states at each of 12 time values for a total of 48 logical qubits. The physical energy gap between the two oscillating states can be extracted as E 2 − E 1 = 2ω.
For an indication of how D-Wave hardware performance has been improving, we compare our findings with Fig. 1 of Ref. [75]. In that work, a similar method was used to observe Rabi oscillations by using the previous two generations of D-Wave annealers, the original 2000Q (released in 2017) and the low-noise 2000Q (released in 2019). The authors of Ref. [75] used those machines to attain five time values and six time values respectively in their Rabi oscillation computations. Our study has used the newer Advantage hardware (released in 2020) to attain 12 time values in Fig. 10.

VII. SUMMARY AND OUTLOOK
In this work, we have used a D-Wave quantum annealer to compute several quantities in a non-Abelian gauge theory. Although the annealer is designed to focus on the optimization problem for an Ising model, we have demonstrated computations in SU(2) lattice gauge theory for eigenvalues, eigenvectors, vacuum expectation values, and real-time evolution. Quantum computing is presently in an era of noisy qubits, but the graphs in Figs. 3-10 demonstrate that a quantum annealer can already produce results that are reasonably precise and accurate, at least on lattices having only a few plaquettes.
For the one-dimensional plaquette lattices studied in this work, classical preprocessing readily produced the explicit Hamiltonian matrices, and then translation and reflection symmetries were used to block diagonalize those matrices. The QAE algorithm [63][64][65][66] allowed each of those blocks to be entered directly into the quantum annealer without the need to construct them from products of quantum gates. The D-Wave hardware readily determined respectable results for the lowest eigenvalue and eigenvector from each of the smaller blocks. To extend this success to larger blocks, we developed an adaptive QAE algorithm which is presented in Appendix C. The D-Wave ocean software [67] offers several methods for tuning hardware performance and enhancing results, but we obtained good quality output without these aids, retaining only the two mandatory adjustable quantities: the penalty parameter λ of QAE and the chain strength of D-Wave. Each of our calculations needed only a portion of the physical Hilbert space, and no unphysical Hilbert space was present.
Real-time oscillations between two excited states were also computed on the quantum annealer. The basic approach was to use Kitaev-Feynman clock states [72,73] and we implemented these by combining the QAE with the TEDVP [74]. Besides demonstrating the ability to access time evolution on a quantum annealer, our QAE+TEDVP computation is also a method for measuring the energy splitting between excited states.
Our choice to begin with SU(2) pure gauge theory on a two-plaquette lattice follows the work of Ref. [6], where the same system was studied on a gate-based quantum computer. Our work goes beyond that starting point in two ways: by extending the number of plaquettes from N plaq = 2 to 4 and 6, and by increasing the gauge field truncation from j max = 1 2 to 1 and 3 2 . For the range of gauge couplings in our work, calculations at these sequences of N plaq and j max values indicate that contributions which would have arisen from still larger N plaq and j max values are negligible. All of our Hamiltonians were straightforward to implement on the quantum annealer because each new matrix can be encoded directly, without any decomposition into gates, but this advantage comes at the expense of requiring many more qubits than the gate-based approach of Ref. [6].
Gate-based quantum computers have the ability to store a quantum state more efficiently than classical computers, and this will be very significant for lattices that are large enough to address the intended goals of nuclear and particle physics phenomenology. Indeed, this statestorage scaling advantage is a key motivator for Ref. [6] and for the entire field of quantum computation in lattice gauge theories. If quantum annealers are to compete with gate-based hardware, then the storage scaling issue must be addressed by quantum annealers as they continue to evolve toward a fully universal form of adiabatic quantum computing. That will be an important challenge for future algorithms on future adiabatic hardware but, independently, near-term quantum annealers might be useful alongside gate-based hardware in complementary rather than competitive ways. For example, some algorithms relevant to lattice gauge theory might achieve a speedup sufficient to supersede classical computers [54][55][56][57], making quantum annealers a useful tool during the next several years of study on intermediatesized lattices. To assess the power of quantum annealing, further experiments are called for [76], and our work is a step along this path. As quantum hardware continues to evolve rapidly, we anticipate that mutually complementary roles could emerge for several quantum and classical hardware platforms within the broad scope of lattice gauge theory research.
where our convention is to use m for the left or bottom end of a link and to use m for the right or top end. The other seven vertices have similar expressions, The product of the eight vertex states defines the state of the entire lattice. Notice that we always list the gauge links A through L in alphabetical order so the calculation will be self-consistent and have the correct Clebsch-Gordan phases. The labeling of gauge links chosen in Fig. 1 (even horizontal, then odd horizontal, then vertical) is not required, but it does maintain a convenient pattern among the four plaquette operators during the derivations. The first plaquette operator is where each sum includes only the two terms s i = ± 1 2 . Notice that the subscripts on U F and U I have been interchanged because going counterclockwise around the plaquette means we are going from the m end to the m end on those two links. The effect of an operator U is [6,14] U s,s |j, m, m = where the sums over M and M contain only a single nonzero term because the Clebsch-Gordan coefficients vanish unless M = m + s and M = m + s . Applying 1 to our initial state gives Applying a final state to that result allows all sums to be performed and the answer simplifies to as given in Eq. (6). Results for ψ final | i |ψ initial with i = 2, 3, 4 can be obtained simply by translation symmetry from the i = 1 result or by explicit calculation.

Appendix B: VACUUM SECTOR MATRICES
As explained in Sec. III, three spatial symmetries can be used to block diagonalize each Hamiltonian matrix. One explicit example was provided in Eq. (8). The largest blocks that arise from three additional physics systems are provided here. Recall that the largest block is always the one containing the vacuum state.
For N plaq = 4 and j max = 1 2 , the vacuum block is with the basis states shown beside their corresponding rows. For N plaq = 6 and j max = 1 2 , the vacuum block is . For N plaq = 2 and j max = 1, the vacuum block is  where B is the number of rows in the matrix and K is the number of logical qubits to be used. The quantities q i are the binary variables that we introduced in Sec. IV. The quantity F was defined in Eq. (14) and is the dataset that gets provided directly to the D-Wave quantum annealer. The iteration can be terminated after any quantum annealing step, for example when subsequent z values are no longer providing an improved (i.e. smaller) groundstate eigenvalue.
A python implementation of the AQAE algorithm is provided in Ref.

ACKNOWLEDGMENTS
We are grateful to D-Wave Systems Inc. (Burnaby, Canada) for providing us with computing time on their quantum annealing hardware. Our work was supported in part by a Discovery Grant from the Natural Sciences and Engineering Research Council (NSERC) of Canada. S.A. received additional funding from a Dean's Under-graduate Research Award, E.M. from a Carswell Graduate Scholarship, and S.P. from an NSERC Undergraduate Student Research Award.