Quantum Gate Sets for Lattice QCD in the strong coupling limit: $N_f=1$

We derive the primitive quantum gate sets to simulate lattice quantum chromodynamics (LQCD) in the strong-coupling limit with one flavor of massless staggered quarks. This theory is of interest for studies at non-zero density as the sign problem can be overcome using Monte Carlo methods. In this work, we use it as a testing ground for quantum simulations. The key point is that no truncation of the bosonic Hilbert space is necessary as the theory is formulated in terms of color-singlet degrees of freedom (``baryons'' and ``mesons''). The baryons become static in the limit of continuous time and decouple, whereas the dynamics of the mesonic theory involves two qubits per lattice site. Lending dynamics also to the ``baryons'' simply requires to use the derived gate set in its controlled version.


Introduction
Recently, there has been considerable interest generated by proposals for quantum simulations of gauge theories in high-energy physics.This interest has spread to lattice gauge theory (LGT), and in particular, to lattice quantum chromodynamics (LQCD).This development promises to allow direct access to real-time physics as well as to the physics of dense nuclear matter from LQCD.Considerable progress has been made in the universal gate-based approach both in the methods needed [1][2][3][4] as well as their applications to theories with both discrete and continuous gauge groups to varying levels of approximation. 1  Although from current estimates it seems that near-term, large-scale simulations of LQCD on quantum computers are infeasible [6], one can still define physically interesting problems which can be studied in the noisy intermediate-scale quantum (NISQ) era [7].In the gate-based approach, one is of course forced to work in the Hamiltonian formulation.In the case of LGT, the Hamiltonian formulation appeared quite early [8] but was of secondary interest due to the success of Markov chain Monte Carlo methods which 1 See [5] for a recent experimental proposal to realize non-Abelian LGTs using analog quantum simulations.
rely on the Lagrangian formulation.Another pioneering study identified the mapping of operators appearing in the Hamiltonian into spin operators for continuous gauge groups assuming a cutoff in the infinite-dimensional Hilbert space of the gauge bosons [9].This line of analysis has also been extended to fermions [10].In order to alleviate the systematic effects of Hilbert space truncation, one can attempt to improve the Kogut-Susskind Hamiltonian by adding terms which reproduce the desired physics for a given cutoff.This has very recently been done for the strong-coupling limit using a similarity renormalization group approach [11].
The strong-coupling limit of lattice QCD has been a long-standing field of interest, both in the Wilson and staggered formulations.Of particular interest are the reformulation of the theory in color-singlet degrees of freedom and their first numerical simulations [12,13].This work was mainly motivated by the possibility to either cure or alleviate the sign problem, and hence study LGT systems at non-zero baryon density.Equipped with new algorithms, these studies were revisited some years later, leading to a complete mapping of the phase diagram for one staggered flavor in the strong coupling limit (β = 0) [14] and beyond [15].All of the above-mentioned studies start from the Lagrangian formulation and thus, in principle, have little to say about real-time evolution of quantum states.Recently, a continuous-time partition function has been formulated leading to the construction of a Hamiltonian describing the physics of the strong-coupling limit with N f = 1 flavors of massless staggered fermions [16].Further efforts were undertaken to extend this to the case of N f = 2 [17,18].
There exists then ample motivation to study the strong-coupling limit of lattice QCD from the perspective of quantum simulations: For one, as the gauge degrees of freedom can be integrated out exactly in any spatial dimension, Gauss' law is fulfilled implicitly.
Secondly, the resulting variables form a discrete set of physical, color-singlet degrees of freedom (mesons and baryons) governed by a local, sparse nearest-neighbor Hamiltonian.
Finally, due to the mildness of the sign problem in the Euclidean formulation, cross checks between the formulations can be performed in the extended parameter space of finite μ, T at continuous Euclidean time.We restrict ourselves to the case of massless quarks which is conceptually simpler [19] and exhibits an exact remnant chiral symmetry.It therefore seems natural to map the strong-coupling Hamiltonian onto the degrees of freedom relevant for quantum simulations, starting with the one-flavor, β = 2N c /g 2 = 0 case which is the subject of our current study.

Background
The Hamiltonian we study is derived from the continuous Euclidean time limit of strongcoupling lattice QCD with staggered fermions.Strong-coupling lattice QCD can be formulated in terms of dual variables in order to tackle the sign problem, with the full μ B -T phase diagram accessible through Monte Carlo simulations.The partition function of the discrete system for one flavor of massless staggered fermions is given by where S F = x γ η 0 (x) χ(x)e a τ μ q U 0 (x)χ(x + 0) -χ (x + 0)e -a τ μ q U † 0 (x)χ(x) is the staggered quark action S F on an anisotropic lattice of extent N τ × L d .As the temperature is given by T = (a(β)N τ ) -1 on an isotropic lattice, it requires anisotropic lattices to continuously vary the temperature at fixed inverse gauge coupling.At strong coupling β = 0, the lattice spacing a is maximally coarse, while the temporal lattice spacing a τ will be sent to zero.The non-perturbative relation between the bare anisotropy γ , introduced above as a prefactor to favor temporal over spatial fermion hops, and the physical anisotropy ξ ≡ a a τ , has been determined non-perturbatively for various N τ and towards the continuous time limit N τ → ∞ [20].The resulting temperature in lattice units is aT = ξ (γ ) N τ , which simplifies in the continuous time limit to aT = κ γ 2 N τ with κ = 0.7971(3).This is close to the mean-field result where κ = 1.
Owing to the limit of infinite gauge coupling, the gauge link integration DU in Eq. ( 1) can be done exactly [13], followed by the Grassmann integration.This results in the partition function of a dimer model [16,21] with the (3 + 1)-dim.lattice volume decomposed into the disjoint union of mesonic sites M and baryonic sites B .The integers k μ ∈ {0, . . ., 3} are dimers, the selfavoiding loops denote baryons worldlines.Note that the sign σ ( ) ∈ {+1, -1} can be negative for nontrivial baryon loops, but σ ( ) = 1 for static baryons.The normalization N (γ ) can be disregarded.The number of spatial meson hops depends on the temperature, but remains finite at fixed aT as N τ → ∞.This results in the suppression of spatial meson hops of orders higher than γ -2 in units of a t , simplifying Eq. (3) significantly: with M and B denoting the spatial mesonic and baryonic sublattices respectively, which no longer depend on τ ∈ {0, . . ., N τ -1}.The mesonic contribution has non-trivial weights whenever, at some time index τ , a meson hop between nearest neighbors x, x + î occurs.
Here v ( x,τ ) depends on the temporal dimers k ± 0 just before/after the spatial dimer location, and due to k - 0 + k + 0 + 1 = 3, the only possible transitions are In between spatial meson exchange, the temporal dimers at a given site x form states that can be identified as meson occupations, m = {0, π, 2π, 3π}, which correspond to a conserved pion current.For the baryonic contribution, spatial hops are suppressed and we have used μ B /T = a τ μ B N τ for the baryon chemical potential μ B = 3μ q .
One finally takes the limit N τ → ∞ of Eq. ( 4), by substituting γ2 with aTN τ , such that the pion exchange between nearest neighbors can be written as an exponential, and the vertices can be expressed with matrices Ĵ± which raise or lower the meson occupation number, resulting in the following partition function: with the Hamilton and number operators x, y The corresponding ladder operators Ĵx = ( Ĵ+ x ) T are: The trace in Eq. ( 5) is over the Hilbert space H h = x∈ H h x , where the local Hilbert space H h x is 6 dimensional, with states |h = |m, b ∈ H h x and a basis |h i of 6 hadronic states per spatial lattice site, The nature and dimensionality of these states can be explained by the fact that the gauge integration in Eq. ( 1), carried out analytically in the limit of infinite gauge coupling (β = 0), results in the reformulation of the theory of Eq. ( 2) in terms of color-singlet degrees of freedom, "baryons" and "mesons".The subsequent Grassmann integration and continuous time limit constrains the site occupation to 0, . . ., N c "pions", N c = 3, alternatively a site may be occupied by a "baryon" or "anti-baryon". 2 By redefining the meson occupation numbers as |s = |m -3/2 , a particle-hole symmetry becomes evident, as the mesonic creation and annihilation operators fulfill the following algebra: This corresponds to a 4-dim.representation of SU(2) with |j, m ≡ |3/2, s .This construction can be extended to a N f = 2 Hamiltonian [17,18], resulting in: where the Hamiltonian is now composed of a sum of 4 contributions, one for each pion In addition to the baryon number operator NB , the N f = 2 partition function also contains an isospin number operator NI .The number of distinct hadronic states |h is now 92, of which 50 states are purely mesonic.It turns out that all the matrix elements in Ĵ± Q i are positive, allowing one to study both non-zero baryon and isospin chemical potential.
Another possible extension of the Hamiltonian Eq. ( 6) is to perform the strong-coupling expansion [22] before taking the continuous time limit, which results in only temporal plaquettes contributing.Both the N f = 2 formulation and the inclusion of the leadingorder gauge corrections to the Hamiltonian will be discussed in a forthcoming publication.

Implementation and results
The Hamiltonian Ĥ and the number operator N , as defined in Eq. ( 6), have a blockdiagonal structure, with the states |h x of the local Hilbert space being a direct sum of mesonic and baryonic occupational states, |h x = |m x ⊕ |b x .As a practical consequence, this allows us to qubitize the Hamiltonian Eq. ( 6) in the mesonic sector and baryonic sector independently.We begin with the former.

Mesonic gate sets
The four-dimensional local meson states |m i ∈ {0, π, 2π, 3π} can be encoded with two qubits.The Hamiltonian Ĥ represents a nearest-neighbor meson hopping term, seen through the action of the raising and lowering operators, Ĵ+ and Ĵ-, respectively, on the local Hilbert space.Such a local, low-dimensional Hamiltonian can be qubitized with standard methods [23] as we summarize below for the case of 1+1-dim.The method, of course, is completely general and can be applied to higher dimensions.
We start by partitioning the lattice Hamiltonian into mutually commuting even and odd parts where Ĥx e/o both have the nearest-neighbor nn-form of Eq. ( 6), and the labels "even" and "odd" refer to the bonds of the one-dimensional lattice.Using the Pauli decomposition H x e/o = i c i P i into Pauli strings P i ∈ {I, X, Y , Z} ⊗4 , the remaining tasks are that of partitioning the P i into sets ("families") of mutually commuting terms, which can then be simultaneously diagonalized.The Pauli decomposition yields eighteen four-qubit terms (two qubits for the neighboring sites x and y, respectively) which can easily be grouped into four families. 3The partition of the operators into four families Fam k takes the form where the big endian convention is used.The individual terms are multiplied by the coefficients c (k) j which are listed in the Appendix (Eq.( 17)).Each of the Fam k can be diagonalized simultaneously using the tableau formalism described in [25], yielding the sets Diag Fam 1/2 = {IIZZ, IZZZ, ZIIZ, ZZIZ}, The transformations V k which diagonalize the families of operators are given below in gate form.The remaining step in the derivation of the primitive gate sets for the H e/o , is the optimization of the exponentiated sum of diagonal terms exp (i j c(k Here, optimization refers to the number of CNOT gates.Making use of the identity Z a ⊗ Z b = CNOT a,b I a ⊗ Z b CNOT a,b and taking into account mutual cancellations between neighboring CNOTs, we arrive at the gate set for each of the four exponentiated families.As an example, the relevant gates for the exponentiation of the first family are displayed in Fig. 1, where the remaining sets are given in the Appendix, Fig. 6.With the observation that the partitioning of the lattice Hamiltonian into mutually commuting sets of terms in nn-form (Eq.( 12)) is easily generalized to higher dimensions, it can readily be seen that the derivation Eqs. ( 13)-( 15), leading to the gate sets Figs. 1 and  6, holds in any dimension.We restrict ourselves to hypercubic lattice geometries which naturally arise in the context of staggered fermions [8].

Inclusion of baryons
Extending our theory from gauge group U(3) to SU(3) requires the extension of |h = |m to |h = |m, b , as defined in Eq. (8).It is important to note that owing to the limit of infinite gauge coupling and continuous time, baryons remain static in our toy theory.Given an initial state |ψ 0 = ⊗ x |h x , the time evolution of bosonic and baryonic sites decouples (Bottom) The diagonal gate diag_famI of the above unitary, decomposed into elementary operations Figure 2 An example of the two-qubit controlled version of the mesonic Trotter gate where the qubits (q 0 q 1 q 2 ) and (q 3 q 4 q 5 ) have been arranged for readability as [H, N ] = 0, with baryons evolving trivially, only leaving an imprint on the surrounding "pion bath" through a fixed excluded volume.From the practical side of quantum simulations this implies that the gate set identified so far is sufficient to capture the dynamics of the theory with gauge group SU(3) in the zero-temperature limit.
At non-zero chemical potential or temperature, μ B , T > 0 we will need to extend our register to three qubits per site x to represent all six classical states |h i ∈ {0, π, 2π, 3π, B + , B -}, with two states being redundant.As Ĥx e/o and Nx act on this Hilbert space through a direct sum, we can use the three qubits in the register (q 0 q 1 q 2 ) and encode the hadronic sector (meson or baryon) of h x in, say, q 0 , thus having it act as control bit.A mesonic Trotter step thus consists of the gates derived in the last section, Fig. 1 and Fig. 6, but in their controlled version (e.g. on "0"), with one control bit per site.The unitary exp (-iδt Ĥx e/o ) acting on pairwise sets of qubits encoding |h x and |h y of nn-pairs x, y , respectively, is hence represented by the symbolic circuit Fig. 2, where the gate depth and CNOT count are given in Table 1.
The baryon evolution is governed by a single operator in diagonal form (i.e.Z-term in the Hamiltonian) locally at a given site x through the term exp (-it x ωx ) and hence does not need to be Trotterized.We thus only have to take into account the local control bit denoted by q 0 (e.g. on "1"), which encodes the hadronic sector.Furthermore, in order to restrict the action of exp (-itω x ) to the (B + , B -) T part of |h x (encoded on q 2 ), we can add Figure 3 Baryonic evolution gate for a single site.The local control bits q 0 and q 1 encode the hadronic sector and redundant state space, respectively a control by q 1 .This excludes the two redundant states of the 2 3 = 8 possible classical states.This is depicted in Fig. 3, where one can see the only baryonic gate used in the time evolution of the full (SU(3)) theory.

Trotterized time evolution in d = 1 + 1
Approximating the time evolution by a first order product formula, the Trotterized timeevolution operator of our theory is as follows with a systematic error of order O(δt 2 ) for every Trotter step δt with δt = t/N .The symbolic circuit corresponding to the staggered layers of mesonic Trotter gates, followed by a baryonic evolution for a lattice of dimension d = 1 + 1 with linear extent L = 4 is depicted in Fig.
a superposition of states with non-zero mesonic and baryonic occupation numbers.The two layers of a single mesonic Trotter step result from the nn-decomposition of our lattice Hamiltonian (Eq.( 12)).Increasing the dimensionality will add two layers per space dimension, corresponding to the increase in coordination number of our hypercubic lattice.The gate depth of a mesonic Trotter step hence grows linearly with increasing lattice dimension and is unaffected by the lattice volume N (cf.Table 1 for a single mesonic gate), whereas the number of qubits grows linearly, n q = 3N .From the point of view of computational complexity this resource requirement can be compared to the qubit count of LGT in the local multiplet basis [2] at lowest truncation order p = 1 where we have n q ≈ 2N log 2 ( p + 1) for pure Yang-Mills theory, i.e. without matter.
To verify the correctness of our gate decomposition summarized in Fig. 4, we show in Fig. 5 a comparison of exact and Trotterized observables as a function of time t, with the initial state |ψ 0 on a lattice of L = 4 sites (n q = 12 qubits) with periodic boundary conditions.The displayed results correspond to a noise-free classical simulation, assuming infinite shot statistics.Using these time-evolved states, we have computed the timedependent expectation values of the mesonic observables Ĵ1/2 (left), defined in Eq. ( 9) as well as the overlap | ψ 0 |ψ(t) | 2 (right) for several values of the baryon chemical potential aμ B .We observe good agreement with the exact results.

Scalability
It is instructive to compare our implementation proposal for quantum simulations of lattice QCD in the strong coupling limit to existing benchmarks at scale, s.a. the recently published study of lattice QED in d = 1 + 1 (i.e. the Schwinger model) with more than 100 qubits [26] on an IBM Heron processor.There, using a second order Trotter formula to implement the time evolution of the system, 14 Trotter steps were reached with total  CNOT-depth of 370.If we restrict ourselves to the case of a purely mesonic theory (thus using the gauge group U(3) instead of SU(3)) and take our resource requirements given in Table 1 at face value, time evolution of our d = 1 + 1 dimensional theory with a layering as given in Fig. 4 will require a CNOT-depth of 88 per first order Trotter step, thus already permitting 4-5 (depending on transpilation) Trotter steps on current hardware with comparable quantum volume as the mesonic theory requires 2 qubits per staggered lattice site.While this may seem less favorable, we emphasize that going to larger spatial dimension in our theory will only result in a linear increase in CNOT-depth without modification of the proposed gate set.
Turning to the full theory (gauge group SU(3)), the pertinent question is the hardware realization of multi-qubit-controlled gates, in view of the 2-qubit controlled version of the mesonic and baryonic gates, Figs. 2 and 3. Commonly, a decomposition is used [27] which for 2-qubit control would essentially triple the CNOT-depth per Trotter step.More promising would be the use of multi-qubit gates [28] on trapped-ion quantum computers which were demonstrated for the case of 3-qubit and 4-qubit gates with promising fidelity [29].
For actual quantum simulations our circuits displayed in Fig. 1 and Fig. 6 clearly will require the use of error mitigation techniques [30] s.a.zero noise extrapolation (ZEN) [31] and probabilistic error cancellation (PEC) [32] to counter the effect of noise in the form of decoherence and limited gate fidelities, in particular those of the 2-qubit entangling gates.To convert the a priori unknown noise channel into Pauli noise channels, Pauli-twirling can be applied [33].An additional hardware limitation, for a given spatial dimension, is the required qubit connectivity introduced by the dimensionality of the Hilbert space requiring 3 qubits per site to be coupled to the qubits of neighboring sites.

Conclusion and outlook
In summary, we have mapped the Hamiltonian for strong-coupling lattice QCD with one flavor of staggered quarks to qubit degrees of freedom by writing down the complete gate set to simulate its dynamics on a quantum computer.This formulation serves as a preparatory step to actual studies of the model on quantum hardware, where it is evident that this will not only allow one to study the dynamics of the model but also its thermal and finite-density properties, using e.g.thermal pure quantum states [34,35] or the variational methods described in [36] and [37].Strong-coupling lattice QCD is studied on classical computers via the Worm algorithm, introduced in [38] and first applied to SU(3) in [14], which is based on the high-temperature expansion of the partition function.It is very efficient at high and intermediate temperatures, but not at low temperatures.Here, quantum simulations might be expected to outperform classical algorithms relying on importance sampling.
Repeating this exercise for N f > 1 is another avenue to pursue.Considering the state multiplicity diagram given in [18] for N f = 2, one sees that the purely bosonic sector of the theory (B = 0) will require the use of six qubits per lattice site and in the case a baryon is present (|B| = 1) will require five qubits.As the mesonic ladder operators Ĵ± remain block-diagonal and commute with the number operator, the mesonic trotter steps will not change the baryon number.While the gate set for the N f = 1 Hamiltonian is well suited for superconducting quantum computers as the scalability analysis has shown, the quantum gate sets for N f = 2 may requires a different strategy: It may be advantageous to map the system to an effective theory or to make use of a qudit as a multi-level computational unit [39], suggesting that a quantum hardware based on trapped ions may be preferable.It should also be mentioned that our approach has limitations: the dual formulation for staggered fermions requires an unimproved fermion action, and as it does not involve a fermion determinant, we cannot implement rooting.
The most physically relevant extension of the current work would be the inclusion of gauge corrections to the Hamiltonian formulation.This would allow one to make contact to other effective theories in the strong-coupling limit [11], and in the longer term, approach the continuum limit a → 0. The dimension of the N f = 1 Hilbert space at finite β becomes 16-dimensional.Higher order corrections will mainly contribute to the matrix elements of the ladder operators, resulting in a larger number of mutually commuting families.While the leading order gauge correction may still be suited for superconducting qubits, higher orders may result in a CNOT depth, for which other computational devices are beneficial.It should be stressed that both extensions of our Hamiltonian formulation still result in finite-dimensional Hilbert spaces, given N f and a fixed order in the strong coupling expansion, with no further truncations required.

Figure 1 (
Figure 1 (Top) Four qubit quantum gate corresponding to the Fam 1 unitary in the mesonic Trotter step.(Bottom) The diagonal gate diag_famI of the above unitary, decomposed into elementary operations

4 .
Here we have used open boundary conditions for illustrative purposes and chose the initial state |ψ 0 = ⊗ L-1 i=0 (|+ |0 |+ ), where at each lattice site one the following state in the hadron occupation basis |h x = 1 2 (

Figure 4 A
Figure4A symbolic circuit for two mesonic Trotter steps (yellow), followed by a baryonic evolution (blue) for the state |ψ 0 = (|+ |0 |+ ) ⊗L with linear extent L = 4 and open boundary conditions.Each of the mesonic Trotter steps corresponds to a sequential application of the gates defined in Fig.1and Fig.6, where the two-qubit control is visible in the first layer

Figure 5 (
Figure 5 (Left) Comparison of the classically simulated Trotterized expectation value of the mesonic observables Ĵ1/2 for Trotter step size δt = 0.4 (diamonds) and the exact results (dashed lines).(Right)The wave function overlap | ψ 0 |ψ(t) | 2 for several values of the baryon chemical potential aμ B for Trotter step size δt = 0.2.We used a system of linear extent L = 4 and periodic boundary conditions for the initial state |ψ 0 = (|+ |0 |+ ) ⊗L

Table 1
Gate depth and CNOT count of a single mesonic Trotter gate acting on a pair of qubit registers, as displayed in Fig.2