# Digital quantum simulation of a (1+1)D SU(2) lattice gauge theory with ion
qudits

Giuseppe Calajó  Istituto Nazionale di Fisica Nucleare (INFN), Sezione di
Padova, I-35131 Padova, Italy  Giuseppe Magnifico  Dipartimento di Fisica e
Astronomia “G. Galilei”, via Marzolo 8, I-35131 Padova, Italy  Dipartimento di
Fisica, Università di Bari, I-70126 Bari, Italy.  Istituto Nazionale di Fisica
Nucleare (INFN), Sezione di Bari, I-70125 Bari, Italy  Claire Edmunds
Universität Innsbruck, Institut für Experimentalphysik, Technikerstraße 25a,
Innsbruck, Austria  Martin Ringbauer  Universität Innsbruck, Institut für
Experimentalphysik, Technikerstraße 25a, Innsbruck, Austria  Simone Montangero
Dipartimento di Fisica e Astronomia “G. Galilei”, via Marzolo 8, I-35131
Padova, Italy  Padua Quantum Technologies Research Center, Universitá degli
Studi di Padova  Istituto Nazionale di Fisica Nucleare (INFN), Sezione di
Padova, I-35131 Padova, Italy  Pietro Silvi  Dipartimento di Fisica e
Astronomia “G. Galilei”, via Marzolo 8, I-35131 Padova, Italy  Padua Quantum
Technologies Research Center, Universitá degli Studi di Padova  Istituto
Nazionale di Fisica Nucleare (INFN), Sezione di Padova, I-35131 Padova, Italy

###### Abstract

We present a quantum simulation strategy for a (1+1)D SU(2) non-abelian
lattice gauge theory, a hardcore-gluon Hamiltonian Yang-Mills, tailored to a
six-level trapped-ion qudit quantum processor, as recently experimentally
realized [1]. We employ a qudit encoding fulfilling gauge invariance, an SU(2)
Gauss’ law. We discuss the experimental feasibility of generalized Mølmer-
Sørensen gates used to efficiently simulate the dynamics. We illustrate how a
shallow circuit with these resources is sufficient to implement scalable
digital quantum simulation of the model. We also numerically show that this
model, albeit simple, can dynamically manifest physically-relevant properties
specific to non-abelian field theories, such as baryon excitations.

Gauge theories are a fundamental theoretical framework in physics playing an
important role in many active areas of research spanning from high energies
[2] to condensed matter [3] and quantum information science [4]. Their
formulation on a lattice, known as lattice gauge theories, (LGT) [5, 6, 7], is
particularly suited to study non-perturbative effects. Monte Carlo techniques
have been extremely successful over the years [8] in tackling various models,
including quantum chromodynamics, but they are limited by the sign problem to
certain physical regimes, and struggle to capture the real-time evolution
outside equilibrium [9]. In the last years, following the advances in quantum
and quantum-inspired computation, alternative new possibilities to face these
problems have emerged based either on tensor network numerical techniques [10,
11, 12] or on analog and digital quantum simulation [13, 14, 15], involving
different experimental platforms such as cold atoms [16, 17, 18, 19, 20, 21,
22], superconducting circuits [23, 24, 25, 26] and trapped ions [27, 28, 29,
30, 31]. While with the numerical approach both abelian [32, 33, 34, 35, 36,
37, 38, 39, 40, 41, 42] and non-abelian [43, 44, 45, 46, 47, 48] lattice gauge
theory models have been theoretically studied, early demonstrations of quantum
simulation of LGT largely focused on abelian theories [28, 30, 49, 50, 51, 52,
53]. First attempts to tackle non-abelian models have been proposed relying
ether on hybrid quantum computation schemes such as variational eigensolvers
[25], while other quantum simulation encoding proposals are fairly limited in
system size for realistic experimental platforms based on qubits [54, 55, 56,
57].

Figure 1: Sketch of the model.–(a) A 1D Yang-Mills SU(2) lattice gauge model
is encoded in a linear chain of trapped ion qudits. The lattice gauge links
(blue ovals) connecting the matter sites (red circles) are truncated to a
5-dimensional Hilbert space. The mapping is then performed by exploiting the
local dressed basis (5), obtained through the decomposition of the gauge links
into a pair of rishons. (b) Pictorial representation of the local dressed
basis given in Eq. (5). The red and green circles in the links represent the
two colored rishons while the ones in the sites represent the two colored
matter fermions.

One potential avenue to perform large-scale quantum computation relies on
qudit-based quantum processors, which have been proposed in several platforms
such as Rydberg arrays [58, 59], photonic circuits [60] and ultra-cold atomic
mixture [61]. An experimental breakthrough along this direction has been the
demonstration of a universal 7-level optical qudit quantum processor
implemented on a chain of trapped ${}^{40}\textrm{Ca}^{+}$ ions [1]. Using
other ion species, the qudit dimension can be further extended as suggested by
recent results with up to 13 levels [62]. These hardware developments have
already stimulated interesting proposals [63, 64, 65] and experiments [66, 67]
for performing simulations of lattice gauge theories exploring this enlarged
Hilbert space. However, a proof-of-principle experimental demonstration of a
scalable quantum simulation of the dynamics of non-abelian LGTs is still
lacking.

Here, we present a compact rishon representation for a truncated Yang-Mills
SU(2) (1+1)D lattice gauge theory reduced to a 6-dimensional local Hilbert
space, embedding the gauge and fermionic degrees of freedom. This
representation allows us to preserve the non-abelian gauge symmetry while
maintaining the range of the gauge-matter interaction to the nearest
neighbors, in contrast with approaches where the adiabatic elimination of the
gauge degree of freedom requires long-range interaction terms [29]. We present
an experimentally feasible proposal on the 40Ca+ trapped-ion implementation
[1, 66]. In particular, we demonstrate how a digital quantum simulation of the
model can be efficiently implemented by making use of generalized Mølmer-
Sørensen gates (MS) realized by simultaneous driving of multiple transitions
[68]. In this manner, we obtain a shallow circuit for each time step enabling
a feasible digital quantum simulation on current devices capable of
distinguishing between meson and baryon production and thus signaling the non-
abelian nature of the model. Finally, we discuss the major challenges for the
experimental realization of the proposal, but we demonstrate that this is
compatible with current technologies.

This paper is structured as follows. In Sec. I, we introduce the (1+1)D SU(2)
Yang-Mills lattice gauge model and derive the truncated qudit Hamiltonian
defined on a local six-dimensional Hilbert space. In Sec. II, we present the
strategy used to encode the model on a trapped-ion qudits quantum processor
and the generalized Mølmer-Sørensen gates used in the quantum digital
simulation. In Sec. III, we show a few paradigmatic examples of the non-
abelian dynamics inherent in the model. In Sec. IV, we present the strategy
used to perform a quantum digital simulation of the model dynamics. In Sec. V,
we discuss the experimental feasibility of the proposal. Finally, in Sec. VI,
we summarize our results and discuss future directions of research.

##  I SU(2) Lattice Yang-Mills

We start by considering a Hamiltonian Yang-Mills lattice gauge model, with
SU(2) color symmetry, in one spatial dimension, focusing on the low energy
regime. The model illustrates a (flavorless) fermionic matter field coupled to
an SU(2) gauge field, as depicted in Fig. 1(a). The matter field, representing
quarks of bare mass $m_{0}$, is described by staggered fermions
$\hat{\psi}_{na}$ [69], with two colors
$a\in\\{{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r},{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}\\}$
(say red and green), living on the lattice sites $n$, and satisfying standard
Dirac anticommutation rules
$\\{\hat{\psi}_{na},\hat{\psi}^{\dagger}_{n^{\prime}b}\\}=\delta_{n,n^{\prime}}\delta_{a,b}$.
Conversely, the non-abelian gauge field lives on the lattice bonds between
sites $n$ and $n+1$. Following the Kogut-Susskind formulation of gauge fields
on a lattice [6], the system Hamiltonian reads:

$\begin{split}H_{0}=&\frac{c\hbar}{2a_{0}}\sum_{n}\sum_{a,b={\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r},{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}}\left[-i\hat{\psi}^{\dagger}_{na}\hat{U}^{ab}_{n,n+1}\hat{\psi}_{n+1,b}+{\rm
H.c.}\right]\\\
&+m_{0}c^{2}\sum_{na}(-1)^{n}\hat{\psi}^{\dagger}_{na}\hat{\psi}_{na}+g_{0}^{2}\frac{c\hbar}{2a_{0}}\sum_{n}\hat{E}^{2}_{n,n+1}\end{split}$
(1)

where $a_{0}$ is the lattice spacing, $\hbar$ the Planck constant and $c$ the
relativistic speed of free massless particles. The first two terms of Eq. (1)
describe the lattice Hamiltonian of the covariant Dirac equation for massive
quarks. It uses the staggered mass term to create two sublattices,
representing each a component of the two-spinor Dirac field, avoiding the
doubling problem [69]. In this sense, quarks are staggered fermions on even
sites, while anti-quarks are staggered fermion holes on odd sites.

The last term is the pure gauge Hamiltonian, which contains only an electric
component $\hat{E}^{2}_{n,n+1}$ since there are no magnetic fields in one
spatial dimension. The dimensionless coupling $g^{2}_{0}$ depends on the quark
color charge $q_{c}$ and on the lattice spacing $a_{0}$: In one spatial
dimension, it scales as
$g^{2}_{0}=({q_{c}^{2}a_{0}^{2}})/({c\hbar\epsilon_{c}})$ with $\epsilon_{c}$
the vacuum color permittivity, assuming the theory is super-renormalizable.
The SU(2)-electric field energy density is captured by the (dimensionless)
quadratic Casimir operator
$\hat{E}^{2}_{n,n+1}=|\hat{\mathbf{L}}_{n,n+1}|^{2}=|\hat{\mathbf{R}}_{n,n+1}|^{2}$,
where the algebra operators $\hat{L}^{(\nu)}_{n,n+1}$ and
$\hat{R}^{(\nu)}_{n,n+1}$ with coordinates $\nu\in\\{x,y,z\\}$ are
respectively the left and right group generators of the gauge transformation
on each link. The gauge-field algebra is defined by the commutation rules

$\displaystyle\left[\hat{L}^{(\nu)}_{n,n+1},\hat{R}^{(\nu^{\prime})}_{n^{\prime},n^{\prime}+1}\right]$
$\displaystyle=0$ (2)
$\displaystyle\left[\hat{L}^{(\nu)}_{n,n+1},\hat{L}^{(\nu^{\prime})}_{n^{\prime},n^{\prime}+1}\right]$
$\displaystyle=i\delta_{nn^{\prime}}\epsilon^{\nu\nu^{\prime}\nu^{\prime\prime}}\hat{L}^{(\nu^{\prime\prime})}_{n,n+1}$
$\displaystyle\left[\hat{R}^{(\nu)}_{n,n+1},\hat{R}^{(\nu^{\prime})}_{n^{\prime},n^{\prime}+1}\right]$
$\displaystyle=i\delta_{nn^{\prime}}\epsilon^{\nu\nu^{\prime}\nu^{\prime\prime}}\hat{R}^{(\nu^{\prime\prime})}_{n,n+1}$
$\displaystyle\left[\hat{L}^{(\nu)}_{n,n+1},\hat{U}^{ab}_{n^{\prime},n^{\prime}+1}\right]$
$\displaystyle=-\delta_{nn^{\prime}}\sum_{c}\frac{\sigma^{(\nu)}_{ac}}{2}\hat{U}^{cb}_{n,n+1}$
$\displaystyle\left[\hat{R}^{(\nu)}_{n,n+1},\hat{U}^{ab}_{n^{\prime},n^{\prime}+1}\right]$
$\displaystyle=\delta_{nn^{\prime}}\sum_{c}\hat{U}^{ac}_{n,n+1}\frac{\sigma^{(\nu)}_{cb}}{2}$

where $\epsilon^{\nu\nu^{\prime}\nu^{\prime\prime}}$ is the Levi-Civita symbol
for SU(2) and $\sigma^{(\nu)}$ are the Pauli matrices.

$A^{(1)}$ |  $A^{(2)}$ |  $B^{(1)}$ |  $\qquad\qquad\quad B^{(2)}$
---|---|---|---
$\begin{pmatrix}0&&&\sqrt{2}&&\\\ &0&1&&&\\\ &1&0&&&1\\\ \sqrt{2}&&&0&\sqrt{2}&\\\ &&&\sqrt{2}&0&\\\ &&1&&&0\\\ \end{pmatrix}$ |  $\begin{pmatrix}0&&&\sqrt{2}i&&\\\ &0&i&&&\\\ &-i&0&&&i\\\ -\sqrt{2}i&&&0&\sqrt{2}i&\\\ &&&-\sqrt{2}i&0&\\\ &&-i&&&0\\\ \end{pmatrix}$ |  $\begin{pmatrix}0&&-\sqrt{2}i&&&\\\ &0&&-i&&\\\ \sqrt{2}i&&0&&-\sqrt{2}i&\\\ &i&&0&&-i\\\ &&\sqrt{2}i&&0&\\\ &&&i&&0\\\ \end{pmatrix}$ |  $\quad\begin{pmatrix}0&&\sqrt{2}&&&\\\ &0&&1&&\\\ \sqrt{2}&&0&&\sqrt{2}&\\\ &1&&0&&1\\\ &&\sqrt{2}&&0&\\\ &&&1&&0\\\ \end{pmatrix}$
Table 1:  Four of the six matrices that define the effective model on qudits. These four are related to the covariant Dirac transport, and they are factors of a nearest-neighbor interaction $A^{(1)}_{n}B^{(1)}_{n+1}+A^{(2)}_{n}B^{(2)}_{n+1}$.  $M$ |  $C$ |  $D^{(L)}$ |  $\qquad\qquad\qquad D^{(R)}$
---|---|---|---
$\begin{pmatrix}0&&&&&\\\ &0&&&&\\\ &&1&&&\\\ &&&1&&\\\ &&&&2&\\\ &&&&&2\\\ \end{pmatrix}$ |  $\begin{pmatrix}0&&&&&\\\ &2&&&&\\\ &&1&&&\\\ &&&1&&\\\ &&&&0&\\\ &&&&&2\\\ \end{pmatrix}$ |  $\begin{pmatrix}+1&&&&&\\\ &-1&&&&\\\ &&+1&&&\\\ &&&-1&&\\\ &&&&+1&\\\ &&&&&-1\\\ \end{pmatrix}$ |  $\quad\;\begin{pmatrix}+1&&&&&\\\ &-1&&&&\\\ &&-1&&&\\\ &&&+1&&\\\ &&&&+1&\\\ &&&&&-1\\\ \end{pmatrix}$
Table 2:  The last two matrices that define the effective model, related to
mass term (before staggerization) and chromoelectric energy density
respectively. The table also shows the matrices of the link fermion parity
selection rule.

To represent these operators in a matrix form it is useful to express them in
the chromoelectric basis of states $|jm_{L}m_{R}\rangle$, where
$j\in\mathbb{N}/2$ labels the spin irreducible representations (irreps),
$m_{R}\in\\{-j..j\\}$ labels a spin shell $j$, and $m_{L}\in\\{-j..j\\}$
labels a state within the spin shell adjoint to $j$.The gauge field algebra
operators in this basis [70]

$\displaystyle\langle
j^{\prime}m^{\prime}_{L}m^{\prime}_{R}|\hat{E}^{2}|jm_{L}m_{R}\rangle$
$\displaystyle=j(j+1)\delta_{j,j^{\prime}}\delta_{m_{L},m_{L}^{\prime}}\delta_{m_{R},m_{R}^{\prime}}$
(3) $\displaystyle\langle
j^{\prime}m^{\prime}_{L}m^{\prime}_{R}|\hat{U}^{ab}|jm_{L}m_{R}\rangle$
$\displaystyle=\left(C^{j,m_{L}}_{j^{\prime},m^{\prime}_{L};\frac{1}{2},a}\right)^{*}C^{j^{\prime},m^{\prime}_{R}}_{j,m_{R};\frac{1}{2},b}$

where $C^{J,M}_{j_{1},m_{1};j_{2},m_{2}}=\langle
j_{1},m_{1};j_{2},m_{2}|J,M\rangle$ are the Clebsh-Gordan coefficients, i.e.
the fusion rules, for SU(2).

The model (1) is designed to be symmetry-invariant under the gauge
transformations generated by
$\hat{\mathbf{G}}_{n}=\left(\hat{\mathbf{R}}_{n-1,n}+\hat{\mathbf{S}}_{n}+\hat{\mathbf{L}}_{n,n+1}\right)$,
where
$\hat{S}^{(\nu)}_{n}=\frac{1}{2}\sum_{a,b}\sigma^{(\nu)}_{ab}\hat{\psi}^{\dagger}_{na}\hat{\psi}_{nb}$
generates the color-rotations for the quarks. Under this observation, the non-
abelian Gauss’ law, which defines the sector of physical states
$|\Psi_{\text{phys}}\rangle$, reads

$\left|\hat{\mathbf{G}}_{n}\right|^{2}|\Psi_{\text{phys}}\rangle=0\;\;\forall
n\qquad\mbox{(Gauss' Law)}$ (4)

corresponding to the absence of a color background.

###  I.1 Hardcore gluon approximation

To digitally quantum simulate Hamiltonian (1) is necessary to truncate the
infinite local gauge Hilbert space to a finite dimension. This can be done at
low energies by employing the Quantum Link Model (QLM) [71] formalism, where
only a finite set of shells $j$ are considered. Basically, a cutoff is chosen
for the $\hat{E}^{2}$ energy term, and the spin shells $j$ above the cutoff
are discarded. Our proposal for digital quantum simulation considers a QLM
cutoff that includes the $j=0$ and $j=\frac{1}{2}$ shells, that is, the
smallest cutoff which allows quarks to form bound states of baryons and
mesons. In practice in our descriptions we keep all the states that are
reachable from the bare vacuum with a single application of $\hat{U}^{ab}$: in
analogy to cold atom physics, we refer to this cutoff strategy as “hardcore
gluons”, and it is a reasonable approximation for strong coupling regimes [46,
47, 72].

Ultimately, this picture describes a 5-dimensional gauge field Hilbert space
at each bond, spanned by the link basis set
$\\{|00\rangle,|{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r}{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r}\rangle,|{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}\rangle,|{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r}\rangle,|{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r}{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}\rangle\\}$.
Within this representation it is possible to decompose each gauge field bond
into a pair of exotic fermion sub-orbitals (rishons). Then, the parallel
transporter can be efficiently decomposed as
$\hat{U}^{ab}_{n,n+1}=\frac{1}{\sqrt{2}}(\zeta^{a}_{n,n+1})_{L}(\zeta^{b\dagger}_{n,n+1})_{R}$
where each of the exotic fermions $\zeta^{a}$ lives on the left (L) and right
(R) rishon of the same link, as sketched in the second row of Fig. 1(a). Each
rishon mode is ultimately 3-dimensional, spanned by the state
$|0\rangle=|j=0,m=0\rangle$ with even fermion parity
$\langle\hat{P}\rangle=+1$, and the states
$|{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r}\rangle=|j=\frac{1}{2},m=+\frac{1}{2}\rangle,$
$|{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}\rangle=|j=\frac{1}{2},m=-\frac{1}{2}\rangle$
with odd fermion parity $\langle\hat{P}\rangle=-1$. In this representation,
the fermion parity on a link becomes an abelian (gauge $\mathbb{Z}_{2}$)
symmetry of the dynamics, which defines the 5-dimensional space.

By contrast, the matter sites are regular Dirac fermions, thus spanned by the
site basis: $|0\rangle$ the fermion vacuum (even fermion parity),
$|{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r}\rangle=\psi^{\dagger}_{{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r}}|0\rangle$
and
$|{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}\rangle=\psi^{\dagger}_{{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}}|0\rangle$
singly-occupied states (odd parity), and finally
$|d\rangle=\psi^{\dagger}_{{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r}}\psi^{\dagger}_{{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}}|0\rangle$
the doubly-occupied state (even). We can now fuse together the (R)n-1,n
rishon, the matter at site $n$, and the (L)n,n+1 rishon into a unique dressed
site. We now enforce the non-Abelian Gauss’ Law of Eq. (4), resulting in an
effective 6-dimensional gauge-invariant basis for the dressed site defined as
the tensor product of the matter field on a lattice site and of the rishon
states on its left and right. This basis, pictorially sketched in Fig. 1(b),
contains only states with even total fermion parity and explicitly reads

$\displaystyle|1\rangle$ $\displaystyle=|0,0,0\rangle$
$\displaystyle|2\rangle$
$\displaystyle=\frac{|{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r},0,{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}\rangle-|{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g},0,{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r}\rangle}{\sqrt{2}}$
(5) $\displaystyle|3\rangle$
$\displaystyle=\frac{|{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g},{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r},0\rangle-|{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r},{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g},0\rangle}{\sqrt{2}}$
$\displaystyle|4\rangle$
$\displaystyle=\frac{|0,{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r},{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}\rangle-|0,{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g},{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r}\rangle}{\sqrt{2}}$
$\displaystyle|5\rangle$ $\displaystyle=|0,d,0\rangle$
$\displaystyle|6\rangle$
$\displaystyle=\frac{|{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r},d,{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}\rangle-|{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g},d,{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r}\rangle}{\sqrt{2}}\,.$

As such, the model preserves fermion parity at each dressed site [73]. In this
dressed site basis, it is relatively straightforward (see App. A) to rewrite
the hardcore-gluon SU(2) Yang-Mills model from (1) into an effective
Hamiltonian

$H=\sum_{n}\left[\hat{A}^{(1)}_{n}\hat{B}^{(1)}_{n+1}+\hat{A}^{(2)}_{n}\hat{B}^{(2)}_{n+1}\right]+m\sum_{j}(-1)^{n}\hat{M}_{n}+g^{2}\sum_{n}\hat{C}_{n}$
(6)

where we set the energy scale based on the hopping term, by rescaling
$H_{0}\to H=\frac{4\sqrt{2}a_{0}}{c\hbar}H_{0}$, to work in natural units,
with a dimensionless Hamiltonian $H$ and dimensionless couplings
$m=\frac{4\sqrt{2}m_{0}a_{0}c}{\hbar}$, and
$g^{2}=\frac{3\sqrt{2}}{4}g_{0}^{2}$. The $6\times 6$ matrix operators
appearing in this expression are reported in Tables 1 and 2.

While the non-Abelian Gauss’ law has been already enforced in the definition
of a dressed basis, we still have to take into account the Link Law which
arose from splitting a gauge field into two rishons. In the effective model
the link law translates into an abelian selection rule

$\left(\hat{D}_{n}^{(L)}\hat{D}_{n+1}^{(R)}-1\right)|\Psi_{\text{phys}}\rangle=0\;\;\forall
n\quad\mbox{(Link Law)}.$ (7)

This symmetry at each link is protected by the Hamiltonian, thus in principle,
it would be sufficient to satisfy the constraint on the initial state of the
dynamics. However, this symmetry may be disrupted by noise and imperfections
in a digital quantum simulator [28, 30]. In Sec. IV.2 we discuss how to
mitigate this error via a postselection procedure.

Another important symmetry to discuss is the conservation of the total baryon
number $\hat{N_{b}}=\frac{1}{2}\sum_{n}(\hat{M}_{n}-1)$. This quantity can be
controlled by appropriately constructing the starting state of the dynamics
and allows the quantum simulator to explore areas of the phase diagram, with
high baryon density, inaccessible to Monte Carlo simulations due to the sign
problem.

##  II Encoding into trapped-ion qudits

The structure of the truncated gauge-preserving Hamiltonian given in Eq. (6),
defined on a dressed local basis of dimension six, suggests a natural
implementation on a qudit-based quantum processor. In the following, we focus
on an implementation using trapped ${}^{40}\textrm{Ca}^{+}$ ions, as presented
in Ref. [1], but the proposed scheme is versatile and applicable also to other
qudit quantum processors. In the trapped-ion experiment of Ref. [1], each
qudit is encoded within the electronic ground state $S_{1/2}$ and the
metastable excited state $D_{5/2},$ as illustrated in Fig. 2. An external
magnetic field splits the $S_{1/2}$ state into two Zeeman sublevels,
$m_{B}=\pm 1/2,$, and the excited state into six Zeeman sublevels ($m_{B}=\pm
5/2,\pm 3/2,\pm 1/2$). This configuration yields a total of eight accessible
qudit levels connected, considering selection rules, by ten allowed electric-
quadrupole transitions ($\Delta m_{B}=0,\pm 1,\pm 2$). A convenient encoding
of the model given in Eq. (6) onto the ion’s levels involves placing the
states $|3\rangle$ and $|4\rangle$ within the $S_{1/2}$ ground state, while
the remaining states reside in the $m_{B}=\pm 3/2$ and $m_{B}=\pm 1/2$ levels
of the $D_{5/2}$ metastable state, as illustrated in Fig. 2. This choice is
motivated by the observation that all the transitions in the model are
interconnected via the two states: $|3\rangle$ and $|4\rangle$. Generic single
qudit operations among these states can be decomposed into at most $d(d-1)/2$
two-level subspace rotations, where $d$ is the qudit dimension. These
rotations have the form
$R(\theta,\phi)=e^{-i(\theta/2)\hat{\sigma}^{s,s^{\prime}}_{\phi}}$, where the
operators $\hat{\sigma}^{s,s^{\prime}}_{\phi}=|s\rangle\langle s^{\prime}|$
connect the qudit states $s$ and $s^{\prime}$, $\theta$ denotes the rotation
angle, and $\phi$ sets the rotation axis. Importantly, in the model defined in
Eq. (6) the mass and gauge Hamiltonians are diagonal, requiring only $5$
elementary rotations, which can be easily implemented with high fidelity.

The hopping part of the Hamiltonian given in Eq. (6) is more demanding. It
necessitates two-qudit operations, which can be broken down into a sequence of
entangling gates operating on pairs of qudit levels:

$\begin{split}\hat{A}^{(1)}&=\hat{\sigma}^{2,3}_{x}+\sqrt{2}\hat{\sigma}^{1,4}_{x}+\sqrt{2}\hat{\sigma}^{4,5}_{x}+\hat{\sigma}^{3,6}_{x}\\\
\hat{A}^{(2)}&=-(\hat{\sigma}^{2,3}_{y}+\sqrt{2}\hat{\sigma}^{1,4}_{y}+\sqrt{2}\hat{\sigma}^{4,5}_{y}+\hat{\sigma}^{3,6}_{y})\\\
\hat{B}^{(1)}&=\hat{\sigma}^{2,4}_{y}+\sqrt{2}\hat{\sigma}^{1,3}_{y}+\sqrt{2}\hat{\sigma}^{3,5}_{y}+\hat{\sigma}^{4,6}_{y}\\\
\hat{B}^{(2)}&=\hat{\sigma}^{2,4}_{x}+\sqrt{2}\hat{\sigma}^{1,3}_{x}+\sqrt{2}\hat{\sigma}^{3,5}_{x}+\hat{\sigma}^{4,6}_{x}\,.\end{split}$
(8)

Note that with the chosen encoding of the model this decomposition involves
only direct transitions between the $S_{1/2}$ and $D_{5/2}$ states. From this
decomposition, it becomes evident that a single next-neighbor ($n,n+1$) block
of the hopping Hamiltonian requires a total of $32$ entangling gates of the
form
$R_{XY}(\varphi)=e^{-i\varphi\hat{\sigma}^{s_{1},s_{2}}_{x}\otimes\hat{\sigma}^{s_{3},s_{4}}_{y}}$.
These rotations can be implemented using Mølmer Sørensen (MS) gates [74].
However, the considerable gate count imposes limitations on the performance of
a quantum simulation of the model. In the following, we will explore how this
scheme can be enhanced by using native qudit gates based on the simultaneous
driving of multiple transitions.

Figure 2: Model encoding into ions qudit.– The six states needed for the model
(6) are encoded within the eight Zeeman sublevels of the $S_{1/2}$ and
$D_{5/2}$ states of ${}^{40}\textrm{Ca}^{+}$ ions. By encoding the states
$|3\rangle$ and $|4\rangle$ into the $S_{1/2}$ manifold all matrix elements
can be implemented via direct transitions. The boxes on the left show the
driving scheme to implement the $\hat{A}^{(k)}$ and $\hat{B}^{(k)}$ matrices,
$k=1,2$, while the boxes on the right show the implementation of the
$\hat{\alpha}^{(k)}$ and $\hat{\beta}^{(k)}$ matrices. The relative Rabi
frequencies needed for each transition are indicated next to the corresponding
arrows with two different colors, red for $\Omega$ and dark red for
$\sqrt{2}\Omega$ .

###  II.1 Native two-qudit gates

The two hopping blocks, denoted as $\hat{A}^{(k)}_{n}\hat{B}^{(k)}_{n+1}$ with
$k=1,2$, can be directly implemented through a generalized Mølmer-Sørensen
scheme [68] that involves the simultaneous driving of all four direct
transitions in each ion, as depicted in Fig. 2. Similar to standard MS gates,
assume that each transition of interest, $\omega_{s,s^{\prime}}$, is driven by
a pair of lasers with frequencies $\omega^{1}_{L}$ and $\omega^{2}_{L}$,
featuring opposite detunings, i.e.,
$\omega^{1}_{L}=\omega_{s,s^{\prime}}+\delta$ and
$\omega^{2}_{L}=\omega_{s,s^{\prime}}-\delta$. To achieve the correct matrix
elements reported in Tab. 1, two out of the four transitions are driven with a
Rabi frequency of $\Omega$, while the other two are driven with
$\sqrt{2}\Omega$, see Fig. 2. Additionally, the phases $\phi_{n}$ associated
with each ion’s driving are chosen to align with the correct phase pattern. We
assume, to work in the Lamb-Dicke regime, $\eta\ll 1$, with $\eta$ being the
Lamb-Dicke parameter, and to separate the blue and red sidebands, coming from
each laser pair, via a rotating wave approximation valid in the regime where
the lasers are tuned close to the motional sideband with frequency $\nu$,
i.e., $|\delta-\nu|\ll\delta$. A hopping block for two target ions $n$ and
$n+1$ can then be implemented via the generalized MS Hamiltonian

$\bar{H}^{(k)}_{\rm
AB}=\hbar\frac{\eta\Omega}{2}\left[\hat{A}^{(k)}_{n}+\hat{B}^{(k)}_{n+1}\right][\hat{a}^{\dagger}e^{i(\nu-\delta)t}+\hat{a}e^{-i(\nu-\delta)t}],$
(9)

where $\hat{a}$ ($\hat{a}^{\dagger}$) is the phonon annihilation (creation)
operator. In the regime of weak driving, $\eta\Omega\ll|\nu-\delta|$, the
phononic bath dispersively mediates interactions among the two ions according
to the effective Hamiltonian

$\bar{H}^{(k)}_{\rm
AB}\simeq\hbar\frac{(\eta\Omega)^{2}}{4(\nu-\delta)}\left[\hat{A}^{(k)}_{n}+\hat{B}^{(k)}_{n+1}\right]^{2},$
(10)

where, in order to minimize the population of the motional mode, we assumed to
let the system evolve for a time $d\bar{t}=2\pi\ell/|\nu-\delta|$, with $\ell$
being a positive integer [74, 68]. To perform the dynamics in natural units,
as in Eq. (6), we re-scale this time by the rate associated with the MS gate,
i.e. $d\bar{t}\rightarrow dt=\pi(\eta\Omega/|\nu-\delta|)^{2}$. The
Hamiltonian (10), once re-scaled, exactly reproduces the desired hopping block
up to unwanted single qudit rotations coming from the terms
$(\hat{A}^{(k)}_{n})^{2}$ and $(\hat{B}^{(k)}_{n+1})^{2}$. Upon aggregating
contributions from all transitions, these rotations simplify to
straightforward diagonal matrices. These terms can then be combined with the
mass and gauge Hamiltonians of Eq. (6), resulting in the following single
qudit Hamiltonian:

$H_{n}=m(-1)^{n}\hat{M}+g^{2}\hat{C}-H_{A^{2}}-H_{B^{2}}$ (11)

where the diagonal terms in natural units read

$\begin{split}&H_{A^{2}}={\rm diag}(2,1,2,4,2,1)\\\ &H_{B^{2}}={\rm
diag}(2,1,4,2,2,1)\,.\end{split}$ (12)

Note that, for $n=1$, besides the gauge and mass terms, only $H_{A^{2}}$
contributes to Hamiltonian (11), while for $n=N$ only the $H_{B^{2}}$ term
does. The Hamiltonian terms given in Eq. (10) and Eq. (11) constitute the
fundamental building block of the digital quantum simulation to be discussed
in Section IV. Importantly, with this procedure involving generalized MS gates
based on the simultaneous driving of four transitions, only 2 entangling
operations (one per each hopping block) are necessary to implement the hopping
between two neighboring sites.

Figure 3: Model phenomenolgy.–(a) Pair production. Particle density ($\rm
a_{1}$) and the difference between double and single particle occupancy
densities ($\rm a_{2}$) as a function of time for different mass-coupling
values as indicated in the figure. (b) Baryon diffusion. ($\rm b_{1}$) Double
occupancy probability distribution as a function of time capturing the
evolution of two baryons initially excited at positions $n=6$ and $n=14$ with
$g^{2}=m=0.5$. To better resolve the dynamics we subtracted the probability
distribution associated to the evolution of the vacuum state. ($\rm b_{2}$)
The interference pattern of the process described in ($\rm b_{1}$) obtained by
subtracting the individual diffusion of each baryon. (c) String dynamics.
Single ($\rm c_{1}$) and double ($\rm c_{2}$) occupancy probability
distributions as a function of time for a quark-antiquark string initially
excited within the sites $n=8-13$ with $g^{2}=m=5$. Here we again subtracted
the vacuum state evolution. ($\rm c_{3}$) Electric field and single and double
occupancy string densities as a function of time for the same scenario as the
two panels above. All the simulations have been carried out for a lattice of
$N=20$ sites solving the dynamics of the model (6) with matrix product states
setting the truncation tolerance to $\rm tol=10^{-7}$ and the maximum bond
dimension to $\rm D_{\rm max}=200$.

###  II.2 Intermediate scheme

The native qudit-gate scheme outlined in the previous section, offers the
advantage of achieving an exceptionally short circuit depth. Nevertheless, the
implementation of this method presents technical challenges, primarily
stemming from the requirement for fine-tuned calibration of all driven
transitions, see Sec. V. Fortunately, this demand can be released by an
intermediate scheme relying on the simultaneous driving of only two distinct
disjoint transitions. The core idea involves decomposing the interaction
matrices as $\hat{A}^{(k)}=\hat{\alpha}^{(k)}_{1}+\hat{\alpha}^{(k)}_{2}$ and
$\hat{B}^{(k)}=\hat{\beta}^{(k)}_{1}+\hat{\beta}^{(k)}_{2}$, where $k=1,2$ and

$\begin{split}\hat{\alpha}^{(1)}_{1}&=\tilde{\sigma}^{2,3}_{x}+\sqrt{2}\tilde{\sigma}^{1,4}_{x}\qquad\qquad\quad\hat{\alpha}^{(2)}_{1}=-(\tilde{\sigma}^{2,3}_{y}+\sqrt{2}\tilde{\sigma}^{1,4}_{y})\\\
\hat{\alpha}^{(1)}_{2}&=\tilde{\sigma}^{3,6}_{x}+\sqrt{2}\tilde{\sigma}^{4,5}_{x}\qquad\qquad\quad\hat{\alpha}^{(2)}_{2}=-(\tilde{\sigma}^{3,6}_{y}+\sqrt{2}\tilde{\sigma}^{4,5}_{y})\\\
\hat{\beta}^{(1)}_{1}&=\tilde{\sigma}^{4,6}_{y}+\sqrt{2}\tilde{\sigma}^{1,3}_{y}\qquad\qquad\quad\hat{\beta}^{(2)}_{1}=\tilde{\sigma}^{4,6}_{x}+\sqrt{2}\tilde{\sigma}^{1,3}_{x}\\\
\hat{\beta}^{(1)}_{2}&=\tilde{\sigma}^{2,4}_{y}+\sqrt{2}\tilde{\sigma}^{3,5}_{y}\qquad\qquad\quad\hat{\beta}^{(2)}_{2}=\tilde{\sigma}^{2,4}_{x}+\sqrt{2}\tilde{\sigma}^{3,5}_{x}\end{split}\,.$
(13)

For this decomposition, in contrast to the earlier approach involving only two
entangling gates, a total of 8 operations are required to reproduce the
hopping between two neighboring sites. The different contributions can be once
again implemented by simultaneously driving the two target transitions with a
pair of bi-chromatic pulses characterized by Rabi frequencies
($\Omega,\sqrt{2}\Omega$), as illustrated in Fig. 2. Assuming the validity of
the same assumptions employed in the preceding section, we derive the
effective Hamiltonian

$\bar{H}^{(k)}_{\rm\alpha_{q}\beta_{q^{\prime}}}\simeq\hbar\frac{(\eta\Omega)^{2}}{4(\nu-\delta)}\left[\hat{\alpha}^{(k)}_{q,n}+\hat{\beta}^{(k)}_{q^{\prime},n+1}\right]^{2},$
(14)

where $q,q^{\prime}=1,2$ and we again set the time step to
$d\bar{t}=2\pi\ell/|\nu-\delta|$. As before, the phonon-mediated interaction
induces unwanted single qudit transitions that reduce, after re-scaling in the
natural units, to twice the same diagonal matrices of in Eq. (LABEL:eq:corr),
$H_{n}=m(-1)^{n}\hat{M}+g^{2}\hat{C}-2H_{A^{2}}-2H_{B^{2}}$.

##  III Model phenomenology

In this section we consider a few paradigmatic examples of non-trivial
dynamics occurring in the model by evolving the initial state
$|\Psi(t=0)\rangle$ under the time-dependent Schrödinger equation ruled by the
Hamiltonian given in Eq. (6). In particular, we identify phenomena and
observables where the non-abelian nature of the theory manifests itself with
distinctive features and that could be implemented incurrent state of the art
experiments.

###  III.1 Vacuum fluctuations

We first consider the phenomena of particle density fluctuations starting from
an initial false vacuum [16, 19, 35]. Let us assume the system to be initially
in the Dirac sea, i.e. the bare vacuum represented by the ground state of the
free part of Hamiltonian (6) for large positive mass and coupling, which
consists in alternated empty and doubly occupied quark and anti-quark sites,
$|\Psi(t=0)\rangle=|5\rangle|1\rangle...|5\rangle|1\rangle$. When the hopping
turns on, this state does not represent the real ground state anymore and it
undergoes a non-trivial dynamical evolution: depending on the Hamiltonian
parameters, spontaneous production of quark anti-quark pairs out of the vacuum
takes place. This effect can be quantified by the lattice particle density
counting for the total number of created quark and anti-quark excitations,

$\rho(t)=\rho_{s}(t)+\rho_{d}(t)\,,$ (15)

where $\rho_{e}(t)=1/(2N)\sum_{n}p_{e}(n,t)$ represent respectively the single
$(e=s)$ and double $(e=d)$ particle occupancy densities computed from the
probability distributions:

$p_{s}(n,t)=|\langle\Psi(t)|3\rangle_{n}|^{2}+|\langle\Psi(t)|4\rangle_{n}|^{2}\,,$
(16)

and

$p_{d}(n,t)=\begin{cases}|\langle\Psi(t)|5\rangle_{n}|^{2}+|\langle\Psi(t)|6\rangle_{n}|^{2}\;\;\;\;\text{if
$n\in\rm even$}\\\
|\langle\Psi(t)|1\rangle_{n}|^{2}+|\langle\Psi(t)|2\rangle_{n}|^{2}\;\;\;\;\text{if
$n\in\rm odd$}\,.\end{cases}$ (17)

In Fig. 3($\rm a_{1}$) we simulated the evolution of the particle density for
different mass-coupling ratios starting from the Dirac sea. For small masses,
$m\ll 1$, particle production is energetically favorable, and after a
transient the system reaches a steady state with approximately one particle
per lattice site. For large masses, $m\gg 1$, instead the particle production
is more expensive and recombination effects let the system oscillate between
the Dirac and the true vacuum. To probe the non-abelian nature of the model we
compute separately the contributions to the particle density coming from the
single and double particle occupancy densities. In turn, these quantities
encode the population of color-singlet pairs formed by a quark and an anti-
quark occupying two neighboring sites with an excited shared gauge link
(meson), and quark-quark (baryon) or anti-quark–anti-quark (anti-baryon) pairs
occupying the same site. Note that the baryon population is absent in abelian
models where just one particle can be hosted in each lattice site. To quantify
the different weights of the two contributions to the total particle density
we plot in Fig. 3($\rm a_{2}$) the difference between double and single
particle occupation densities as a function of time for the same mass-coupling
ratios as in Fig. 3($\rm a_{1}$). The figure shows how single particle
occupation is dominant in the weak coupling regime, $g^{2}\ll 1,m$, while in
the strong coupling regime $g^{2}\gg 1,m$ the production of baryons is more
energetically favorable. Such baryon production has no analogs in abelian
theories and is a clear signature of the underlying $SU(2)$ nature of the
model.

Figure 4: Digital quantum simulation.–(a)-(b) Left panels: Circuit
decomposition for the two strategies outlined in the main text, (a) full and
(b) intermediate scheme. In both $\mathcal{D}$ indicates the circuit depth.
Right panels: Dynamics of the particle density evolving from the Dirac vacuum
for $N=3$ lattice sites obtained through the model (6) (continuous lines) and
via a first-order Suzuki-Trotter evolution of the illustrated circuits. In the
latter, we explicitly included the phononic mode as in Eq. (9) truncating at
the eighth level. The corresponding fidelity between the exact and the
digitally simulated state is also shown. The time steps used for the Suzuki-
Trotter evolution are indicated in the figures. In all the simulations, we set
$m=g^{2}$.

###  III.2 Baryon diffusion

The second dynamical phenomenon under consideration involves the diffusion of
two pairs of quarks (baryons) initially localized on two even sites of the
lattice, with the rest of the system being in the Dirac vacuum. In the strong
coupling limit, where $g\gg m,1$, these pairs play the role of composite
particles in the lattice, hopping from one matter site to another, weakly
exciting all other allowable configurations. This regime, characterized by a
fourth-order process with a rate of $J_{\rm eff}=16/g^{4}m$, occurs on a
timescale significantly longer than the natural unit (see App.B).
Consequently, it is not practically feasible for digital quantum simulation
due to the requirement for numerous time steps (see Sec.IV).

To investigate baryon diffusion on a shorter timescale, we examine the non-
perturbative regime where $g\sim m\lesssim 1$. Specifically, we focus on the
scenario where two baryons are initially excited out of the Dirac vacuum,
$\hat{S}^{n_{1}}_{51}\hat{S}^{n_{2}}_{51}|{\rm GS}\rangle$, where
$\hat{S}^{n}_{ss^{\prime}}=|s\rangle\langle s^{\prime}|_{n}$ with $n_{1}$ and
$n_{2}$ being even numbers. In Fig. 3($\rm b_{1}$), the evolution of the
double occupancy probability distribution, $p_{d}$, is plotted, with the
vacuum fluctuation subtracted to reveal the diffusion process. The two baryons
disperse within a wavefront cone with an aperture of approximately $\pm 8t$,
governed by the quark hopping rate. The interference of the two baryons is
depicted in Fig. 3($\rm b_{2}$) and is obtained by subtracting the individual
baryon propagations.

###  III.3 String dynamics

Finally, we consider the dynamics of an initially excited string. In this
context, string-breaking is a typical phenomenon occurring in confined LGTs
where an initial string excitation undergoes evolution by breaking into pairs
of particles and antiparticles. This phenomenon has been extensively studied
numerically in (1+1)D U(1) models [33, 35] and in a (1+1)D SU(2) theory using
a different rishon representation than the one employed here [43]. To
initialize the state configuration, we apply the following string operator of
length $l$ to the Dirac sea vacuum:
$\mathcal{\hat{S}}=\hat{\psi}^{\dagger}_{n_{s}}\hat{U}_{n_{s},n_{s}+1},...\hat{U}_{n_{s}+l-1,n_{s}+l}\hat{\psi}_{n_{s}+l}$
where we assume $n_{s}$ to be even, corresponding to the creation of a quark,
and $n_{s}+l$ to be odd, corresponding to the creation of an anti-quark on
that site. Expressed in the local dressed basis of Eq. (5) this string reads:

$\mathcal{\hat{S}}=\hat{S}^{n_{s}}_{41}\hat{S}^{n_{s}+1}_{65}\hat{S}^{n_{s}+2}_{21}\;...\;\hat{S}^{n_{s}+l-2}_{65}\hat{S}^{n_{s}+l-1}_{21}\hat{S}^{n_{s}+l}_{35}\,.$
(18)

Here we focus again on distinctive non-abelian features by distinguishing
between meson and baryon resonant production from the string. To observe
resonant production of baryons and mesons, we consider the case where the
initial string of energy $E_{\rm str}=2m+2lg^{2}$ resonates with the energy of
$l+1$ baryons, each having energy $E_{\rm bar}=2m$, and $(l+1)/2$ mesons, each
having energy $E_{\rm mes}=2m+2g^{2}$. This resonance condition is fulfilled
for both processes when $g^{2}=m$. We then let the system evolve in time with
fixed mass and coupling set to large values to ensure that most of the energy
remains confined within the string, preventing quick dispersion toward the
system edges. In Fig. 3($\rm c$) we indeed observe resonant excitations in
time of the single and double occupancy probability distributions within the
string, signaling the creation of the string of baryons and mesons. These
resonant oscillations are better resolved in Fig. 3($\rm c_{3}$) where we plot
the single $(e=s)$ and double $(e=d)$ occupancy densities averaged within the
string, $\tilde{\rho}_{e}(t)=1/(2(l+1))\sum^{n_{s}+l}_{n=n_{s}}p_{e}(n,t)$,
along with the electric field string density defined as

$\begin{split}&\tilde{\rho}_{E}=\frac{1}{2l}\sum^{n_{s}+l-1}_{n=n_{s}}|\langle\Psi(t)|2\rangle_{n}|^{2}+|\langle\Psi(t)|4\rangle_{n}|^{2}+|\langle\Psi(t)|6\rangle_{n}|^{2}\\\
&+|\langle\Psi(t)|2\rangle_{n+1}|^{2}+|\langle\Psi(t)|3\rangle_{n+1}|^{2}+|\langle\Psi(t)|6\rangle_{n+1}|^{2}\,,\end{split}$
(19)

which accounts for the number of rishons on each link.

##  IV Digital quantum simulation

In this section, we explore how the dynamics previously described can be
simulated digitally using a trapped-ion qudit quantum processor. To digitally
simulate the dynamics, we use the first-order Suzuki-Trotter decomposition:

$U(t)=e^{-iHt_{\rm f}}\simeq\left(\prod_{j}e^{-ih_{j}dt}\right)^{n_{\rm
ST}}\,,$ (20)

where, $dt=t_{\rm f}/n_{\rm ST}$ denotes the time step ruling the precision of
the Suzuki-Trotter evolution, $n_{\rm ST}$ is the number of Suzuki-Trotter
steps, and $t_{\rm f}$ is the final time of the evolution. The terms $h_{j}$
represent the Hamiltonian building blocks composing Hamiltonian (6), which we
presented in Sec. II. The evolution of the model is then obtained by
executing, for each time step, the unit circuit cell (in parallel across the
lattice), depicted in Fig. 4 for both schemes presented in Sec. II. The unit
circuit cell involves three lattice sites and is composed of the generalized
MS gates, implementing the hopping part of the model, followed by a single-
qudit rotation on each ion encompassing the gauge and mass terms of the
Hamiltonian (11), and the correcting rotations. Considering that single-qudit
rotations can be performed with high fidelity, we define the circuit depth of
the model, $\mathcal{D}$, as the number of generalized MS gates per unit
circuit cell. Then the total number of two-qudit gates, $N_{\rm gates}$,
required to perform the entire simulation on a chain of $N$ sites scales
linearly with the system size and the Trotter steps:

$N_{\rm gates}=\frac{(N-1)}{2}n_{\rm ST}\mathcal{D},$ (21)

with $\mathcal{D}=4$ and $\mathcal{D}=16$ respectively for the first and
second scheme presented in Sec. II. Note that the total gate count of our
proposal for both schemes is drastically reduced with respect to a computation
based on a decomposition into two-level entangling gates (see Sec. II), which
has $\mathcal{D}=64$.

Figure 5: Reproducing model phenomenology.–(a) Difference between double and
single particle occupancy densities as a function of time for a lattice with
$N=3$ sites, considering various mass-coupling values as depicted in the
figure. (b) Baryon diffusion is analyzed for a lattice with $N=3$ sites and
$g^{2}=m=0.5$. The probabilities of double occupancy for the first and third
atoms, denoted as $p_{d}(1,t)$ and $p_{d}(3,t)$, along with the single
occupancy probability for the second atom, denoted as $p_{s}(2,t)$, are
studied over time. The initial condition involves an anti-baryon excited in
the first site. (c) Electric field, as well as single and double occupancy
string densities, are plotted over time for a lattice comprising $N=4$ sites
with $g^{2}=m=5$. The string is initially excited within the sites $n=1-4$. In
all plots, continuous lines represent the exact evolution of the model (6),
while discrete points correspond to the first-order Suzuki-Trotter evolution,
with time step $dt=0.01\pi$, performed using (9), which includes 8 levels of
the phononic mode.

For a realistic numerical simulation of the model, we consider small lattice
sizes, specifically $N=3$ and $N=4$. These sizes are chosen as they retain the
essential characteristics of the SU(2) dynamics under consideration while
being accessible with a short circuit implementation. The generalized MS gates
are simulated explicitly including the vibrational degree of freedom using
Hamiltonian (9) and its analogue for the intermediate scheme.

###  IV.1 Reproducing the system dynamics

To assess the performance of the proposed quantum digital simulation, we
compare the results of the exact expected dynamics with those derived from the
digital scheme. We apply this comparison to the same illustrative examples of
dynamics discussed in Section III. As a first example, we consider in Fig. 4
the particle density evolution of the Dirac vacuum on $N=3$ lattice sites for
the two schemes outlined in Sec. II. To perform a comparison with the exact
evolution, we also compute the state fidelity,
$\mathcal{F}(t)=\langle\Psi(t)|\varrho_{\rm q}(t)|\Psi(t)\rangle$, where
$\varrho_{\rm q}(t)$ represents the reduced density operator of the qudit
system after tracing over the phononic degrees of freedom, and
$|\Psi(t)\rangle$ is the state evolved under the Schrödinger equation. The two
schemes exhibit comparable performances with the precision of the simulation
ruled by the Suzuki-Trotter step $dt$. Importantly, even for large time steps,
$dt=0.04\pi$, the initial density peak, arising from particle production and
characterizing the response time of the system, can be accurately approximated
with just $n_{\rm ST}=3$ Suzuki-Trotter steps. Higher simulation precision can
be reached by performing a second-order Suzuki-Trotter decomposition, as
presented in App. C and discussed in Sec. V.

Despite the limited lattice size, the digital simulation effectively captures
distinct production rates for paired (doubly occupied sites) and unpaired
(singly occupied sites) particles. This is exemplified in Fig. 5(a) where we
observed a similar phenomenology as the one depicted in Fig. 3(a) for the
disparity between double and single particle occupation densities. In addition
to particle production from the vacuum, the digital simulation, when confined
to a few lattice sites, successfully replicates the other two phenomena
discussed in Section III: baryon diffusion and string dynamics. The former is
illustrated in Fig. 5(b) where we considered the hopping of an anti-baryon
excited on the first site towards the third site. This process is quantified
by observing the double particle occupancy probability on the first lattice
site, involving an intermediate single particle population on the second site.
To enhance the resolution of this effect, we subtracted the Dirac vacuum
evolution from the simulation, consistent with the approach taken in Section
III. The final phenomenon under consideration is string breaking. For this
scenario, we extend our analysis to a slightly larger system size,
specifically $N=4$. Proceeding as in Fig. 3($\rm c_{3}$), we calculate the
electric field, as well as the densities of single and double occupancy
strings, for a string initially spanning the entire lattice. As in Fig. 3($\rm
c_{3}$) we set the string energy to be in resonance with meson and baryon
production. This tuning leads to pronounced oscillations in the quark and
electric field densities, serving as indicators of the resonant excitations of
mesons and baryons.

###  IV.2 Link parity preservation

As discussed in Sec. I, the employed rishon representation introduces an extra
$\mathbb{Z}_{2}$ symmetry, which we set to an even number of rishons per link.
While this symmetry is maintained throughout the Hamiltonian evolution, it may
be compromised in actual simulations due to experimental errors leading to
leakage from this subspace. These errors can be compensated in post-selection
by rejecting states that violate the symmetry, leading to a potentially
significant reduction in simulation errors.

To select which data to retain we make use of the parity matrices $D^{(L)}$
and $D^{(R)}$ defined in Table 2, which ensures the link fermion parity
selection rule. When applied jointly to two neighboring sites, these operators
yield a positive outcome if the link maintains the correct rishon parity and a
negative outcome if the link deviates from the correct parity sector.
Specifically, $\hat{D}^{(L)}_{n}\hat{D}^{(R)}_{n+1}|\Psi_{\rm
phys}\rangle=|\Psi_{\rm phys}\rangle$ and
$\hat{D}^{(L)}_{n}\hat{D}^{(R)}_{n+1}|\Psi_{\rm unphys}\rangle=-|\Psi_{\rm
unphys}\rangle$, where $|\Psi_{\rm phys}\rangle$ and $|\Psi_{\rm
unphys}\rangle$ represent the physical and unphysical states, respectively.
These operators establish a truth table that facilitates the post-selection of
measurements conforming to the rishon parity.

We propose two different schemes for post-selecting the data, based
respectively on destructive and non-destructive measurements. The first
exploits the truth table established by the parity operators to determine
which states to retain at the end of the simulation. This selection relies
solely on the analysis of the qudit populations at each site, rejecting
configurations involving adjacent unphysical states–those with just a single
rishon on the link. However, it’s important to note that this procedure does
not entirely eliminate the possibility that, in previous time steps, the
system may have exited the correct parity sector and returned to it afterward.
Nevertheless, it serves as a filter and, since it only requires the qudit
populations for each time step during the readout, it does not incur
additional computational costs. The second protocol exploits an ancilla qubit
to conduct state-preserving measurements. This approach offers the advantage
of continuously monitoring the rishon parity during the time evolution but
demands higher resources. The details of this scheme are discussed in App. D.

##  V Experimental considerations

In this section, we explore the experimental challenges associated with
realizing a quantum digital simulation. Several sources of error can impact
the accuracy and effectiveness of the simulation. These include deviations
from the Lamb-Dicke regime and the rotating-wave approximation, motional
heating, and magnetic field fluctuations [68, 75, 76]. Another significant
hurdle is the implementation of the generalized MS gates via simultaneously
driving multiple transitions. Achieving the correct target operation requires
the precise calibration of many (linear in the number of transitions) coupled
control parameters. Most notably, optical Stark shifts are induced by each
laser tone [77, 78], which creates an additional layer of complexity. In the
following, we discuss some of these experimental challenges, which must be
overcome for the successful experimental realization of the proposed quantum
digital simulation.

###  V.1 Magnetic field fluctuations

To evaluate the impact of external magnetic field fluctuations, we assumed
infinitely correlated noise in time and space: throughout each single
realization of the dynamics, the magnetic field $\mathcal{B}^{z}$ perceived by
the ion qudits is constant and uniform. Therefore, the produced Zeeman shifts
correspond to $H_{Z0}=-\mathcal{B}^{z}\sum_{n}\hat{\mathbf{m}}^{z}_{n}$ with
the magnetic dipole operator

$\mathbf{m}^{z}=\frac{1}{2}\begin{pmatrix}-3\mu_{D}&&&&&\\\ &-\mu_{D}&&&&\\\
&&-\mu_{S}&&&\\\ &&&\mu_{S}&&\\\ &&&&\mu_{D}&\\\ &&&&&3\mu_{D}\\\
\end{pmatrix},$ (22)

where we consider a ratio $w=\mu_{D}/\mu_{S}\approx 0.6$ between the magnetic
moments of the D5/2 and the S1/2 spin-orbitals, compatible with Ref. [1]. Once
we express this Hamiltonian shift in natural units, just as we did for the
model Hamiltonian in Eq. (6), we can similarly identify a magnetic field
realization via the dimensionless parameter
$b=\frac{2\sqrt{2}a_{0}\mu_{s}}{c\hbar}\mathcal{B}^{z}$. Then, the rescaled
(dimensionless) Hamiltonian $H_{Z}=\frac{4\sqrt{2}a_{0}}{c\hbar}H_{Z0}$ simply
reads $H_{Z}=-b\sum_{n}\hat{F}_{n}$ with

$F={\rm diag}(-3w,-w,-1,1,w,3w)\,,$ (23)

and we incorporated the magnetic fluctuations Hamiltonian into the Suzuki
Trotter evolution of the Dirac vacuum depicted in Fig. 4(a). We averaged over
100 realizations of time evolution, while randomly choosing in each
realization $r$ the (static, uniform) rescaled magnetic coupling $b_{r}$ from
a uniform (flat) distribution within the interval $b_{r}\in[-\Delta b,\Delta
b]$. For a comparative analysis with the exact evolution, we computed the
state fidelity, $\mathcal{F}(\tilde{t})$, at the time $\tilde{t}$,
corresponding to the initial peak in the particle density observed in Fig.
4(a). This fidelity is shown in Fig. 6(a) as a function of the strength of the
magnetic field fluctuations, $\Delta b$, for both first and second-order (see
App. C) Suzuki-Trotter evolutions, considering the same Trotter steps as in
Fig. 4. The figure demonstrates that the fidelity is not significantly
impacted for magnetic fluctuations inducing level shifts up to $10\%$ of the
particle hopping rate, set by $(\eta\Omega)^{2}/(2|\nu-\delta|)$ in the
digital simulation. For expected hopping rates on the order of $0.1-1$ kHz,
this implies a requirement for magnetic field fluctuations to be maintained
below the threshold of nT. Achieving this level of precision can be realized
with current magnetic shielding techniques [79].

###  V.2 Impact of other errors in the simulation

To estimate the impact of other sources of error in the digital simulation we
define the Fidelity of a single generalized two-qudit MS gate as
$\mathcal{F_{\rm MS}}$. The overall performance of the simulation is then
given by

$\mathcal{P}(t)=\mathcal{F}(t)[\mathcal{F_{\rm MS}}]^{N_{\rm gates}}$ (24)

where the total number of gates $N_{\rm gates}$ is defined in Eq. (21). In
Fig. 6(b) we plot this simulation performance as a function of different
values of the two-qudit gate fidelity and to the number of Suzuki-Trotter
steps necessary to reach, as before, the peak of the particle density creation
at time $\tilde{t}$. In particular, we compare the results obtained with a
first and a second-order Suzuki-Trotter evolution, with the latter having
$\mathcal{D}=6$ circuit depth (see App. C). This estimation shows that, with a
two-qudit gate fidelity of the order of $\mathcal{F_{\rm MS}}\gtrsim 0.99$ it
is possible to obtain a good simulation performance, $\mathcal{P}\sim 90\%$,
sufficient for probing the dynamics of the model. Such two-qudit gate
fidelities are achievable under the assumption that the errors associated with
each transition are independent of each other. With this assumption, the
overall error grows linearly with the number of transitions with respect to
the standard MS qubit gate infidelity, which, according to the state of the
art, can reach values of the order of $0.2\%$ [80]. However, in the actual
implementation, the generalized MS gate performance will depend on the
accumulated errors coming from each driven transition, which generically are
mutually dependent. The final gate fidelity then will rely on an optimal
calibration process, as discussed in the next section.

Figure 6: Errors estimation.–(a) State infidelity calculated with respect to
the expected state at time $\tilde{t}=0.377$ and the digitally simulated one
as a function of the strength of the magnetic field fluctuations for the same
configuration as in Fig. 4(a). The results were obtained using the first
scheme given in Fig. 4(a) and the continuous and dashed lines represent
respectively the results obtained with a first (ST1) and second order (ST2)
Suzuki Trotter evolution. The grey lines indicate the ideal performance
obtained in the absence of a magnetic field. We consider the same Suzuki-
Trotter time steps as in 4, which require respectively $n_{\rm ST}=12$ and
$n_{\rm ST}=3$ steps to reach $\tilde{t}$. (b) Estimation of the simulation
performance to reach the peak in the particle density at $\tilde{t}=0.377$
versus the number of Trotter steps and the fidelity of a single generalized MS
gate. The left and right panels are obtained respectively for the first and
second-order Suzuki-Trotter evolution.

###  V.3 Calibration challenges

Using multiple driving fields, even a standard qubit MS gate requires the
calibration of four control parameters. In contrast to local gates, these
control parameters are non-linearly correlated and can thus not be calibrated
independently. As an example, consider a bichromatic light field,
symmetrically detuned from an optical transition, as described in Sec. II.
Changing the amplitude of one of the two laser tones non-linearly shifts the
transition frequency, which in turn changes the effective detuning and thus
coupling strength of the other tone. Moreover, due to the multi-level
structure of the ion, these Stark shifts are different for each of the optical
transitions. For a qubit MS gate, calibration typically involves four control
parameters that can be calibrated iteratively or in an automated fashion using
Bayesian techniques [81].

For generalized MS gates, as discussed in Sec II, the situation becomes
significantly more difficult. Not only does a change in one parameter affect
the entangling dynamics between the two states connected by the addressed
transition, but it now also affects the dynamics of states coupled by a
simultaneously driven transition. While manual calibration of such a coupled
multi-parameter landscape seems very challenging, Bayesian parameter
estimation techniques might be extendable to this scenario, given a suitable
cost function. Notably, the second approach discussed in Sec. II mediates
these challenges, since it requires simultaneous driving only on disjoint
transitions. This makes it possible to reduce the coupling between the
parameters and pre-calibrate the operations on the two transitions
independently to a large extent.

Finally, we observe that a further difficulty for the calibration process
comes from the fact that the two-qudit gates presented in Section II require
to drive two different sets of transitions on the two ions $n$ and $n+1$. This
issue can be solved by applying a rotation of the local dressed basis via
unitary transformations, as discussed in Appendix E.

##  VI Conclusions

In summary, we have considered a Yang-Mills SU(2) 1D lattice gauge theory with
dynamically coupled matter truncated at the lowest levels exhibiting non-
trivial dynamics and non-abelian features. We introduced a compact rishon
representation, embedding local gauge and fermionic degrees of freedom within
a dimension-six Hilbert space, which has a natural encoding onto a qudit
quantum processor based on 40Ca+ trapped ions. We presented three different
schemes relying on entangling gates based respectively on a single, double
disjoint and four simultaneously driven transitions per ion. Relevantly, for
the latter, we demonstrated that an efficient digital quantum simulation of
the model can be accomplished with a notably short circuit depth. This result
is facilitated by harnessing the computational advantages inherent in qudits
and is based on the efficacy of generalized Mølmer-Sørensen gates. Note that,
as also pointed out in other works [63, 64, 65, 66], just the simple encoding
of the model into qudits already brings substantial advantages with respect to
qubit-based quantum digital simulations, which usually require larger
computation resources and circuit depth and the capability of engineering
long-range and/or three or more body interactions [25, 82, 55, 57, 83, 84,
85].

The hardcore gluon approximation can in principle be relaxed while keeping an
analogous encoding [46], effectively softening the truncation of the gauge
fields. However, such an operation would require the inclusion of additional
electronic levels in the programmable atomic quantum processing, potentially
achievable with different atomic species [62]. This proposal provides an
experimentally feasible and potentially scalable pathway for observing non-
abelian lattice gauge phenomena, such as those from high-energy physics or
emerging in condensed matter models, in currently available qudit quantum
processors.

##  VII Acknowledgments

The authors thank Guido Pagano, Michael Meth, and Jad Halimeh for valuable
discussions. They acknowledge financial support from: the European Union via
QuantERA2017 project QuantHEP, via QuantERA2021 project T-NiSQ, via the
Quantum Technology Flagship project PASQuanS2, and via the NextGenerationEU
project CN00000013 - Italian Research Center on HPC, Big Data and Quantum
Computing (ICSC); the Italian Ministry of University and Research (MUR) via
PRIN2022-PNRR project TANQU, and via Progetti Dipartimenti di Eccellenza
project Frontiere Quantistiche (FQ); the WCRI-Quantum Computing and Simulation
Center (QCSC) of Padova University. The authors also acknowledge computational
resources by Cloud Veneto and the ITensors Library for the MPS calculations
[86, 87]. The UIBK team acknowledges funding by the European Union under the
Horizon Europe Programme–Grant Agreement 101080086– NeQST, by the European
Research Council (ERC, QUDITS, 101080086), and by the Austrian Science Fund
(FWF) through the SFB BeyondC (FWF Project No. F7109) and the EU-QUANTERA
project TNiSQ (FWF Project No. N-6001). Views and opinions expressed are
however those of the author(s) only and do not necessarily reflect those of
the European Union or European Climate, Infrastructure and Environment
Executive Agency (CINEA). Neither the European Union nor the granting
authority can be held responsible for them.

##  Appendix A Rishon representation

Once the hardcore gluon approximation is adopted, the parallel transporter
$U^{ab}$ reads, in the chromoelectric 5-dimensional basis of the gauge field
space [70],

$\hat{U}^{ab}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{@{\hspace{0.1ex}}c@{\hspace{0.1ex}}|@{\hspace{0.1ex}}c@{\hspace{0.2ex}}c@{\hspace{0.2ex}}c<EMAIL_ADDRESS>0.43057pt\lx@intercol\hfil 0\hfil\hskip 0.43057pt\vrule\hskip
0.43057pt&+\delta_{a{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r}}\delta_{b{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}}\hfil\hskip
0.86108pt&-\delta_{a{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r}}\delta_{b{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r}}\hfil\hskip
0.86108pt&+\delta_{a{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}}\delta_{b{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}}\hfil\hskip
0.86108pt&-\delta_{a{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}}\delta_{b{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r}}\hfil\hskip
0.43057pt\\\ \hline\cr\hskip
0.43057pt\lx@intercol\hfil-\delta_{a{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}}\delta_{b{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r}}\hfil\hskip
0.43057pt\vrule\hskip 0.43057pt&0\hfil\hskip 0.86108pt&\hfil\hskip
0.86108pt&\hfil\hskip 0.86108pt&\hfil\hskip 0.43057pt\\\ \hskip
0.43057pt\lx@intercol\hfil-\delta_{a{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}}\delta_{b{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}}\hfil\hskip
0.43057pt\vrule\hskip 0.43057pt&\hfil\hskip 0.86108pt&0\hfil\hskip
0.86108pt&\hfil\hskip 0.86108pt&\hfil\hskip 0.43057pt\\\ \hskip
0.43057pt\lx@intercol\hfil+\delta_{a{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r}}\delta_{b{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r}}\hfil\hskip
0.43057pt\vrule\hskip 0.43057pt&\hfil\hskip 0.86108pt&\hfil\hskip
0.86108pt&0\hfil\hskip 0.86108pt&\hfil\hskip 0.43057pt\\\ \hskip
0.43057pt\lx@intercol\hfil+\delta_{a{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r}}\delta_{b{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}}\hfil\hskip
0.43057pt\vrule\hskip 0.43057pt&\hfil\hskip 0.86108pt&\hfil\hskip
0.86108pt&\hfil\hskip 0.86108pt&0\hfil\hskip 0.43057pt\\\ \end{array}\right).$
(25)

Here we briefly decompose it as a pair of exotic fermion operators, each one
acting on a sub-orbital (rishon). First of all, we need a practical strategy
to define exotic fermion operators: A valid fermionic operator has a local
action $\tilde{F}$ that inverts a local parity $P=P^{\dagger}=P^{-1}$, so that
$\\{P,\tilde{F}\\}=0$. In analogy to the Jordan-Wigner transformation, the
global action of the fermion operator $\hat{F}$ is the string

$\hat{F}_{j}=\dots P_{j-2}\otimes
P_{j-1}\otimes\left(\tilde{F}\right)_{j}\otimes\mathbb{1}_{j+1}\otimes\mathbb{1}_{j+2}\dots,$
(26)

in contrast to a boson $\hat{B}$ (or spin) global action which reads instead

$\hat{B}_{j}=\dots\mathbb{1}_{j-2}\otimes\mathbb{1}_{j-1}\otimes\left(\tilde{B}\right)_{j}\otimes\mathbb{1}_{j+1}\otimes\mathbb{1}_{j+2}\dots.$
(27)

This prescription ensures that any two fermion operators on different modes
anticommute as they should $\\{\hat{F}_{j},\hat{F}_{j^{\prime}\neq
j}\\}=\\{\hat{F}_{j},\hat{\psi}^{(\dagger)}_{j^{\prime}\neq j}\\}=0$. In this
formalism, it is straightforward to define a Dirac fermion lattice field

as well as a Majorana one

where the subscript ‘$F$’ specifies that the operator is meant to be
understood as a fermion, i.e. that his global action has a parity string
attached as in Eq. (26).

We can then define a pair of exotic fermion operators for the rishon space,
written in the basis
$\\{|0\rangle,|{\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r}\rangle,|{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}\rangle\\}$,
namely

with the corresponding local fermion parity

$P_{\zeta}=\quantity(\begin{array}[]{c|cc}1&0&0\\\ \hline\cr 0&-1&0\\\
0&0&-1\\\ \end{array}).$ (42)

One can then check that the parallel transporter can be decomposed as

$\hat{U}^{ab}_{j,j+1}=\frac{1}{\sqrt{2}}\left(\zeta_{a}\right)_{jL}\left(\zeta_{b}^{\dagger}\right)_{j+1,R}.$
(46)

$Q_{L}$ |  $\qquad\qquad\qquad Q_{R}$
---|---
$\begin{pmatrix}0&&&\sqrt{2}&&\\\ &0&1&&&\\\ &&0&&&1\\\ &&&0&\sqrt{2}&\\\ &&&&0&\\\ &&&&&0\\\ \end{pmatrix}$ |  $\quad\;-i\begin{pmatrix}0&&\sqrt{2}&&&\\\ &0&&1&&\\\ &&0&&\sqrt{2}&\\\ &&&0&&1\\\ &&&&0&\\\ &&&&&0\\\ \end{pmatrix}$
Table 3:  The effective model matrices $\hat{Q}_{L}$ and $\hat{Q}_{R}$ written
in the gauge-invairant canonical basis of the dressed site.

Performing such decomposition greatly simplifies the covariant hopping from
the lattice Yang-Mills Hamiltonian, specifically

$\displaystyle H_{\text{hop}}$
$\displaystyle=\frac{c\hbar}{2a_{0}}\sum_{j}\sum_{a,b={\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r},{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}}\left[-i\hat{\psi}^{\dagger}_{ja}\hat{U}^{ab}_{j,j+1}\hat{\psi}_{j+1,b}+{\rm
H.c.}\right]$ (49)
$\displaystyle=\frac{\sqrt{2}c\hbar}{4a_{0}}\sum_{j}\sum_{a,b={\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r},{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}}\left[-i\hat{\psi}^{\dagger}_{ja}\hat{\zeta}_{jLa}\hat{\zeta}^{\dagger}_{j+1,Rb}\hat{\psi}_{j+1,b}+{\rm
H.c.}\right]$
$\displaystyle=\frac{\sqrt{2}c\hbar}{4a_{0}}\sum_{j}\left[\hat{Q}_{Lj}^{\dagger}\hat{Q}_{R,j+1}+\hat{Q}_{Lj}\hat{Q}_{R,j+1}^{\dagger}\right]$

where the operators
$\hat{Q}_{Lj}=\sum_{a={\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r},{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}}\hat{\zeta}^{\dagger}_{jLa}\hat{\psi}_{ja}$
and
$\hat{Q}_{Rj}=-i\sum_{a={\color[rgb]{0.8359375,0.2421875,0.30859375}\definecolor[named]{pgfstrokecolor}{rgb}{0.8359375,0.2421875,0.30859375}r},{\color[rgb]{0,0.5,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.5,0}g}}\hat{\zeta}^{\dagger}_{jLa}\hat{\psi}_{ja}$
are explicitly gauge invariant and genuinely local (preserves the fermion
parity on the dressed site). On the 6-dimensional logical basis for the
dressed site, these matrices read as reported in Table 3. They are not
hermitian, however, so they still need to be manipulated for quantum
simulation with trapped ions. Therefore we substitute

$\displaystyle\hat{A}^{(1)}$ $\displaystyle=\hat{Q}_{L}+\hat{Q}^{\dagger}_{L}$
$\displaystyle\hat{B}^{(1)}$ $\displaystyle=\hat{Q}_{R}+\hat{Q}^{\dagger}_{R}$
(50) $\displaystyle\hat{A}^{(2)}$
$\displaystyle=i(\hat{Q}_{L}-\hat{Q}^{\dagger}_{L})\quad$
$\displaystyle\hat{B}^{(2)}$
$\displaystyle=i(\hat{Q}_{R}-\hat{Q}^{\dagger}_{R})$

where we defined four hermitian operators. Now, by noticing that
$\hat{A}^{(1)}\otimes\hat{B}^{(1)}+\hat{A}^{(2)}\otimes\hat{B}^{(2)}$ is equal
to
$2\hat{Q}^{\dagger}_{L}\otimes\hat{Q}_{R}+2\hat{Q}_{L}\otimes\hat{Q}_{R}^{\dagger}$
we can conclude that

$H_{\text{hop}}=\frac{\sqrt{2}c\hbar}{8a_{0}}\sum_{j}\left[\hat{A}_{j}^{(1)}\hat{B}_{j+1}^{(1)}+\hat{A}_{j}^{(2)}\hat{B}_{j+1}^{(2)}\right].$
(51)

Which is the expression for the nearest-neighbor term which appears in the
main text.

Regarding the chromoelectric energy density, we can split the energy
contribution half-half to each rishon (the quadratic Casimir eigenvalue must
be the same on the two halves due to the link law). Basically

$H_{\text{elec}}=g^{2}\frac{c\hbar}{4a_{0}}\sum_{j}\left(|\hat{\mathbf{R}}_{j,j+1}|^{2}+|\hat{\mathbf{L}}_{j,j+1}|^{2}\right),$
(52)

but the quadratic Casimir on the rishon space simply reads

$|\hat{\mathbf{L}}|^{2}\backslash|\hat{\mathbf{R}}|^{2}=\quantity(\begin{array}[]{c|cc}0&0&0\\\
\hline\cr 0&3/4&0\\\ 0&0&3/4\\\ \end{array})=\frac{3}{4}\hat{K}.$ (53)

Therefore, in conclusion

$H_{\text{elec}}=g^{2}\frac{3c\hbar}{16a_{0}}\sum_{j}\left(\hat{K}_{Lj}+\hat{K}_{Rj}\right),$
(57)

and the sum $\hat{K}_{L}+\hat{K}_{R}$ on a dressed site is exactly the
operator $\hat{C}$.

##  Appendix B Perturbative model in the strong coupling limit

In the strong coupling limit, $g\gg 1$, the excitation of the link becomes
energetically costly and the field can be adiabatically eliminated. The
model’s description is then reduced to the bare matter states $|1\rangle$ and
$|5\rangle$, which we rename $|\uparrow\rangle=|5\rangle$ and
$|\downarrow\rangle=|1\rangle$. A simple hopping process between these two
states can be obtained in perturbation theory and leads to the following
Heisenberg model in a staggered magnetic field

$H_{\rm
Heis}=-\frac{4}{\bar{g}^{2}}\sum_{n}\bm{\sigma}^{n}\cdot\bm{\sigma}^{n+1}+m\sum_{n}(-1)^{n}\hat{\sigma}_{z}^{n}$
(58)

where $\bar{g}^{2}=g^{2}+m$ and
$\bm{\sigma}=(\hat{\sigma}_{x},\hat{\sigma}_{y},\hat{\sigma}_{z})$. Within
this model the Dirac vaccuum correspond to the antiferromagnetic ground state.
Assuming to initially flip a single spin out of this vacuum, corresponding to
the excitation of a barion in the original model, a second order hopping
process occur in the limit of $\bar{g}^{2}m\gg 1$. This process gives rise to
the following tight-binding Hamiltonian:

$H_{\rm eff}=J_{\rm
eff}\sum_{n}\left(\hat{\sigma}_{+}^{n}\hat{\sigma}_{-}^{n+1}+\hat{\sigma}_{-}^{n}\hat{\sigma}_{+}^{n+1}\right)$
(59)

with hopping rate $J_{\rm eff}=16/\bar{g}^{4}m$, which describes the diffusion
of barions at a speed $v_{g}=2J_{\rm eff}$. This result was obtained in a two
step perturbation process. A more precise effective description can be
obtained performing directly a $4$th order perturbation theory in the limit
$g,m\gg J$, considering explicitly the two hopping paths accessible to an
initially excited barion:

$\begin{split}&|5\rangle|5\rangle|1\rangle\quad
E=0\qquad\qquad\qquad|5\rangle|5\rangle|1\rangle\quad E=0\\\
&|5\rangle|4\rangle|3\rangle\quad
E=2m+2g^{2}\qquad\;\;\;|5\rangle|4\rangle|3\rangle\quad E=2m+2g^{2}\\\
&|5\rangle|1\rangle|5\rangle\quad
E=4m\qquad\qquad\quad\;|4\rangle|6\rangle|3\rangle\quad E=4g^{2}\\\
&|4\rangle|3\rangle|5\rangle\quad
E=2m+2g^{2}\qquad\;\;\;|4\rangle|3\rangle|5\rangle\quad E=2m+2g^{2}\\\
&|1\rangle|5\rangle|5\rangle\quad
E=0\qquad\qquad\qquad|1\rangle|5\rangle|5\rangle\quad E=0\\\ \end{split},$
(60)

where we indicated for each state the bare energy cost. Performing an
adiabatic elimination at the fourth order we finally obtain the following
effective coupling rate

$J_{\rm eff}=\frac{16}{m\bar{g}^{4}}+\frac{8}{g^{2}\bar{g}^{4}},$ (61)

which reduces to the one previously derived in the limit of $g\gg m$.

Figure 7: Second order Suzuki-Trotter decomposition.–(a) Unit cell circuit for
the second orders Suzuki-Trotter decompostion. The first and last gate can be
combined in one of duration $dt$ leading to a circuit depth of
$\mathcal{D}=6$. (b) Particle density as function of time for the Dirac vacuum
with $N=3$ lattice sites. The continuous line has been obtained with the exact
model while the dots via the second order Suzuki-Trotter (ST2) evolution of
the illustrated circuit using the time steps indicated in the figure. Below is
shown the corresponding fidelity between the exact and the digitally simulated
state comparing the results of the first (dashed line) and second order
(continuous line) Suzuki-Trotter decomposition. For the ST2 we we included the
phononic mode truncating at the eighth level.

##  Appendix C Second order Suzuki-Trotter decomposition

In the main text we showed how employing a second order Suzuki-Trotter
decomposition can be beneficial for the simulation, even if it implies a
slightly larger circuit depth. The unit circuit cells is presented in Fig.
7(a). The first and last gate can be merged in one during the evolution thus
the the circuit depth is $\mathcal{D}=6$ To ensure in this scheme a reduced
population of the phonic mode we set $d\bar{t}=4\pi/|\nu-\delta|$ for all the
applied gates. In Fig. 7(a) we show the performance of this scheme considering
again the particle creation out of the Dirac vacuum and the corresponding
state Fidelity in time comparing the results obtained with the first and
second order Suzuki-Trotter decomposition. The plot shows a clear improvement
of the simulation performance due to the reduced Suzuki-Trotter error.

##  Appendix D State and link symmetry preserving protocol

The link parity preservation protocol presented in the main text to preserve
the link symmetry has the disadvantage of relying on destructive measurements.
As an alternative scheme, here we propose a strategy based on state-preserving
measurements. This approach offers the advantage of continuously monitoring
(and quantum Zeno protecting) rishon parity during the evolution of the
system. The core concept involves the utilization of an ancilla qubit, which
may also be represented by an unused qudit level, to conduct state-preserving
measurements. Focusing on two adjacent sites, denoted as $n$ and $n+1$, let us
consider a scenario where, during the system’s evolution, the state tends to
leak out of the parity sector due to experimental errors. The generic form of
the state then reads:

$|\psi\rangle=\sum_{n}a_{n}|\psi_{\rm phys}\rangle_{n}+\sum_{n}b_{n}|\psi_{\rm
unphys}\rangle_{n}$ (62)

where $|\psi_{\rm phys}\rangle_{n}$ and $\psi_{\rm unphys}\rangle_{n}$
represent states living in the even and odd symmetry sector. To non-
destructively measure parity, we employ a hybrid version of the ${\rm
CEX}_{c}$ gate, already implemented in the trapped ion qudit quantum processor
[1]

${\rm
CEX}_{c}=\left\\{\begin{array}[]{cc}|c,\uparrow\rangle\leftrightarrow|c,\downarrow\rangle&\\\
|j,k\rangle\rightarrow|j,k\rangle&{\rm if}\;\;j\neq c\end{array}\right.$ (63)

This gate flips the ancilla qubit’s state if the qudit is in the state
$|c\rangle$ and retains the qubit’s state if the qudit is in any other state.
We initialize the ancilla in the state $|\downarrow\rangle_{a}$ and apply a
sequence of ${\rm CEX}_{c}$ gates on each qudit. The choice of the $|c\rangle$
qudit’ s state corresponds to the negative values of the two parity operators
given in Table (2) associated to each qudit, as illustrated in Fig. 8. This
process results in the following state:

$|\psi\rangle=\sum_{n}a_{n}|\psi_{\rm
phys}\rangle_{n}|\downarrow\rangle_{a}+\sum_{n}b_{n}|\psi_{\rm
unphys}\rangle_{n}|\uparrow\rangle_{a}\,.$ (66)

Figure 8: Parity check circuit to maintain the evolution in the correct parity
sector. The details are given in text.

With this procedure, we can monitor the state’s evolution at each time step.
In particular, if the ancilla qubit remains in its initial state, we can be
confident that the evolution has occurred within the correct symmetry sector.
However, if the ancilla qubit flips, indicating a deviation, it suggests a
leakage into the odd symmetry sector. In response to this deviation, we may
applying an $X$ qudit gate defined as:

$X=\sum_{j}|j+1({\rm mod}\;6)\rangle\langle j|$ (67)

on the first qudit. While this operation doesn’t correct errors, it
effectively project the unphysical state back into the correct symmetry
sector. Although not error-correcting per se, this action can assist the
system in maintaining the intended dynamics.

##  Appendix E Local dressed basis rotation

The two-qudit gates discussed in Sec. II require to drive two different set of
transitions on the two ion’s qudit. This is not a fundamental problem per se
but it can affect the complexity of the pulse calibration process. A solution
to this issue consists in performing an unitary transformation of the local
dressed basis on just the even lattice sites using the following unitary
operator:

$\hat{U}=\begin{pmatrix}i&&&&&\\\ &i&&&&\\\ &&0&1&&\\\ &&1&0&&\\\ &&&&-i&\\\
&&&&&-i\\\ \end{pmatrix}.$ (68)

This matrix transforms the operators of the model (6) according to

$\begin{split}&\hat{U}\hat{A}^{(k)}\hat{U}^{\dagger}=-\hat{B}^{(k)}\\\
&\hat{U}\hat{B}^{(k)}\hat{U}^{\dagger}=\hat{A}^{(k)}\\\
&\hat{U}\hat{M}\hat{U}^{\dagger}=\hat{M}\\\
&\hat{U}\hat{C}\hat{U}^{\dagger}=\hat{C}\end{split}$ (69)

with $k=1,2$. The Hamiltonian given in Eq. (6) than reads

$\begin{split}H=&\sum_{n\in\rm
odd}\left[\hat{A}^{(1)}_{n}\hat{A}^{(1)}_{n+1}+\hat{A}^{(2)}_{n}\hat{A}^{(2)}_{n+1}\right]\\\
-&\sum_{n\in\rm
even}\left[\hat{B}^{(1)}_{n}\hat{B}^{(1)}_{n+1}+\hat{B}^{(2)}_{n}\hat{B}^{(2)}_{n+1}\right]\\\
+&m\sum_{j}(-1)^{n}\hat{M}_{n}+g^{2}\sum_{n}\hat{C}_{n}\,.\\\ \\\ \end{split}$
(70)

In this way, as required, the two-qudit gates necessary to implement the
hopping terms always involve the same operator acting on the pair of qudit.
This change of local basis should be applied before and after the sequence of
four two-qudit gates presented in Fig. 4 of the main text. To conclude the
protocol, after having returned to the original basis, the single qudit
compensation operations (see Eq. LABEL:eq:corr) can be applied as in the main
text.

For the disjoint double-transition scheme is also possible to apply local
transformations in order two have two-qudit gates involving the same operator
on both qudit. Compared to the full scheme, the same basis transformation can
not be applied for all the gates but different single qudit rotations should
be applied on each site before and after each two-qudit gate. In particular,
it is possible to reduce all the two qudit gates to be always of the same form
by using the unitary transformations

$\hat{U}_{\alpha^{(k)}_{q}\alpha^{(k^{\prime})}_{q^{\prime}}}\hat{\alpha}^{(k)}_{q}\hat{U}_{\alpha^{(k)}_{q}\alpha^{(k^{\prime})}_{q^{\prime}}}^{\dagger}=\alpha^{(k^{\prime})}_{q^{\prime}}$
(71)
$\hat{U}_{\beta^{(k)}_{q}\alpha^{(k^{\prime})}_{q^{\prime}}}\hat{\beta}^{(k)}_{q}\hat{U}_{\beta^{(k)}_{q}\alpha^{(k^{\prime})}_{q^{\prime}}}^{\dagger}=\alpha^{(k^{\prime})}_{q^{\prime}},$
(72)

with $k,q,k^{\prime},q^{\prime}=1,2$. Again, the single qudit compensation
matrices could be applied similar as in the main text at the end of the series
of gates when the local basis is transformed back to the original one.

## References

  * Ringbauer  _et al._ [2022] M. Ringbauer, M. Meth, L. Postler, R. Stricker, R. Blatt, P. Schindler, and T. Monz, A universal qudit quantum processor with trapped ions, Nature Physics 18, 1053 (2022).
  * Kronfeld [2012] A. S. Kronfeld, Twenty-first century lattice gauge theory: Results from the quantum chromodynamics lagrangian, Annual Review of Nuclear and Particle Science 62, 265 (2012).
  * Fradkin [2013] E. Fradkin,  _Field Theories of Condensed Matter Physics_, 2nd ed. (Cambridge University Press, 2013).
  * Zeng  _et al._ [2019] B. Zeng, X. Chen, D.-L. Zhou, X.-G. Wen, _et al._ ,  _Quantum information meets quantum matter_ (Springer, 2019).
  * Wilson [1974] K. G. Wilson, Confinement of quarks, Phys. Rev. D 10, 2445 (1974).
  * Kogut and Susskind [1975] J. Kogut and L. Susskind, Hamiltonian formulation of wilson’s lattice gauge theories, Phys. Rev. D 11, 395 (1975).
  * Banks  _et al._ [1976] T. Banks, L. Susskind, and J. Kogut, Strong-coupling calculations of lattice gauge theories: (1 + 1)-dimensional exercises, Phys. Rev. D 13, 1043 (1976).
  * Gattringer and Lang [2009] C. Gattringer and C. Lang,  _Quantum chromodynamics on the lattice: an introductory presentation_, Vol. 788 (Springer Science & Business Media, 2009).
  * Calzetta and Hu [2008] E. A. Calzetta and B.-L. B. Hu,  _Nonequilibrium Quantum Field Theory_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2008).
  * Silvi  _et al._ [2014] P. Silvi, E. Rico, T. Calarco, and S. Montangero, Lattice gauge tensor networks, New Journal of Physics 16, 103015 (2014).
  * Dalmonte and Montangero [2016] M. Dalmonte and S. Montangero, Lattice gauge theory simulations in the quantum information era, Contemporary Physics 57, 388 (2016).
  * Silvi  _et al._ [2019a] P. Silvi, F. Tschirsich, M. Gerster, J. Jünemann, D. Jaschke, M. Rizzi, and S. Montangero, The Tensor Networks Anthology: Simulation techniques for many-body quantum lattice systems, SciPost Phys. Lect. Notes , 8 (2019a).
  * Wiese [2013] U.-J. Wiese, Ultracold quantum gases and lattice systems: quantum simulation of lattice gauge theories, Annalen der Physik 525, 777 (2013).
  * Banuls  _et al._ [2020] M. C. Banuls, R. Blatt, J. Catani, A. Celi, J. I. Cirac, M. Dalmonte, L. Fallani, K. Jansen, M. Lewenstein, S. Montangero, _et al._ , Simulating lattice gauge theories within quantum technologies, The European physical journal D 74, 1 (2020).
  * Di Meglio and et al [2023] A. Di Meglio and et al, Quantum computing for high-energy physics: State of the art and challenges. summary of the qc4hep working group (2023), arXiv:2307.03236 .
  * Banerjee  _et al._ [2012] D. Banerjee, M. Dalmonte, M. Müller, E. Rico, P. Stebler, U.-J. Wiese, and P. Zoller, Atomic quantum simulation of dynamical gauge fields coupled to fermionic matter: From string breaking to evolution after a quench, Phys. Rev. Lett. 109, 175302 (2012).
  * Zohar  _et al._ [2015] E. Zohar, J. I. Cirac, and B. Reznik, Quantum simulations of lattice gauge theories using ultracold atoms in optical lattices, Reports on Progress in Physics 79, 014401 (2015).
  * Zohar  _et al._ [2017] E. Zohar, A. Farace, B. Reznik, and J. I. Cirac, Digital lattice gauge theories, Phys. Rev. A 95, 023604 (2017).
  * Kühn  _et al._ [2014] S. Kühn, J. I. Cirac, and M.-C. Bañuls, Quantum simulation of the schwinger model: A study of feasibility, Phys. Rev. A 90, 042305 (2014).
  * Tagliacozzo  _et al._ [2013a] L. Tagliacozzo, A. Celi, A. Zamora, and M. Lewenstein, Optical abelian lattice gauge theories, Annals of Physics 330, 160 (2013a).
  * Tagliacozzo  _et al._ [2013b] L. Tagliacozzo, A. Celi, P. Orland, M. Mitchell, and M. Lewenstein, Simulation of non-abelian gauge theories with optical lattices, Nature communications 4, 2615 (2013b).
  * Yang  _et al._ [2020] B. Yang, H. Sun, R. Ott, H.-Y. Wang, T. V. Zache, J. C. Halimeh, Z.-S. Yuan, P. Hauke, and J.-W. Pan, Observation of gauge invariance in a 71-site bose–hubbard quantum simulator, Nature 587, 392 (2020).
  * Marcos  _et al._ [2013] D. Marcos, P. Rabl, E. Rico, and P. Zoller, Superconducting circuits for quantum simulation of dynamical gauge fields, Phys. Rev. Lett. 111, 110504 (2013).
  * Mezzacapo  _et al._ [2015] A. Mezzacapo, E. Rico, C. Sabín, I. L. Egusquiza, L. Lamata, and E. Solano, Non-abelian su(2) lattice gauge theories in superconducting circuits, Phys. Rev. Lett. 115, 240502 (2015).
  * Atas  _et al._ [2021] Y. Y. Atas, J. Zhang, R. Lewis, A. Jahanpour, J. F. Haase, and C. A. Muschik, Su (2) hadrons on a quantum computer via a variational approach, Nature communications 12, 6499 (2021).
  * Farrell  _et al._ [2024] R. C. Farrell, M. Illa, A. N. Ciavarella, and M. J. Savage, Quantum simulations of hadron dynamics in the schwinger model using 112 qubits, Preprint at 2401.08044 (2024).
  * Hauke  _et al._ [2013] P. Hauke, D. Marcos, M. Dalmonte, and P. Zoller, Quantum simulation of a lattice schwinger model in a chain of trapped ions, Phys. Rev. X 3, 041018 (2013).
  * Martinez  _et al._ [2016] E. A. Martinez, C. A. Muschik, P. Schindler, D. Nigg, A. Erhard, M. Heyl, P. Hauke, M. Dalmonte, T. Monz, P. Zoller, _et al._ , Real-time dynamics of lattice gauge theories with a few-qubit quantum computer, Nature 534, 516 (2016).
  * Muschik  _et al._ [2017] C. Muschik, M. Heyl, E. Martinez, T. Monz, P. Schindler, B. Vogell, M. Dalmonte, P. Hauke, R. Blatt, and P. Zoller, U (1) wilson lattice gauge theories in digital quantum simulators, New Journal of Physics 19, 103020 (2017).
  * Nguyen  _et al._ [2022] N. H. Nguyen, M. C. Tran, Y. Zhu, A. M. Green, C. H. Alderete, Z. Davoudi, and N. M. Linke, Digital quantum simulation of the schwinger model and symmetry protection with trapped ions, PRX Quantum 3, 020324 (2022).
  * Davoudi  _et al._ [2020] Z. Davoudi, M. Hafezi, C. Monroe, G. Pagano, A. Seif, and A. Shaw, Towards analog quantum simulations of lattice gauge theories with trapped ions, Phys. Rev. Res. 2, 023015 (2020).
  * Rico  _et al._ [2014] E. Rico, T. Pichler, M. Dalmonte, P. Zoller, and S. Montangero, Tensor networks for lattice gauge theories and atomic quantum simulation, Phys. Rev. Lett. 112, 201601 (2014).
  * Pichler  _et al._ [2016] T. Pichler, M. Dalmonte, E. Rico, P. Zoller, and S. Montangero, Real-time dynamics in u(1) lattice gauge theories with tensor networks, Phys. Rev. X 6, 011023 (2016).
  * Ercolessi  _et al._ [2018] E. Ercolessi, P. Facchi, G. Magnifico, S. Pascazio, and F. V. Pepe, Phase transitions in ${Z}_{n}$ gauge models: Towards quantum simulations of the schwinger-weyl qed, Phys. Rev. D 98, 074503 (2018).
  * Magnifico  _et al._ [2020] G. Magnifico, M. Dalmonte, P. Facchi, S. Pascazio, F. V. Pepe, and E. Ercolessi, Real Time Dynamics and Confinement in the $\mathbb{Z}_{n}$ Schwinger-Weyl lattice model for 1+1 QED, Quantum 4, 281 (2020).
  * Rigobello  _et al._ [2021] M. Rigobello, S. Notarnicola, G. Magnifico, and S. Montangero, Entanglement generation in $(1+1)\mathrm{D}$ qed scattering processes, Phys. Rev. D 104, 114501 (2021).
  * Emonts  _et al._ [2023] P. Emonts, A. Kelman, U. Borla, S. Moroz, S. Gazit, and E. Zohar, Finding the ground state of a lattice gauge theory with fermionic tensor networks: A $2+1\mathrm{D}$ ${\mathbb{z}}_{2}$ demonstration, Phys. Rev. D 107, 014505 (2023).
  * Funcke  _et al._ [2023] L. Funcke, K. Jansen, and S. Kühn, Exploring the $cp$-violating dashen phase in the schwinger model with tensor networks, Phys. Rev. D 108, 014504 (2023).
  * Angelides  _et al._ [2023a] T. Angelides, L. Funcke, K. Jansen, and S. Kühn, Computing the mass shift of wilson and staggered fermions in the lattice schwinger model with matrix product states, Phys. Rev. D 108, 014516 (2023a).
  * Belyansky  _et al._ [2023] R. Belyansky, S. Whitsitt, N. Mueller, A. Fahimniya, E. R. Bennewitz, Z. Davoudi, and A. V. Gorshkov, High-energy collision of quarks and hadrons in the schwinger model: From tensor networks to circuit qed (2023), arXiv:2307.02522 [quant-ph] .
  * Kebric̆  _et al._ [2023] M. Kebric̆, J. C. Halimeh, U. Schollwöck, and F. Grusdt, Confinement in 1+1d $\mathbb{Z}_{2}$ lattice gauge theories at finite temperature (2023), arXiv:2308.08592 [cond-mat.quant-gas] .
  * Florio  _et al._ [2023] A. Florio, A. Weichselbaum, S. Valgushev, and R. D. Pisarski, Mass gaps of a $\mathbb{Z}_{3}$ gauge theory with three fermion flavors in 1 + 1 dimensions (2023), arXiv:2310.18312 .
  * Kühn  _et al._ [2015] S. Kühn, E. Zohar, J. I. Cirac, and M. C. Bañuls, Non-abelian string breaking phenomena with matrix product states, Journal of High Energy Physics 2015, 1 (2015).
  * Silvi  _et al._ [2017] P. Silvi, E. Rico, M. Dalmonte, F. Tschirsich, and S. Montangero, Finite-density phase diagram of a $(1+1)-d$ non-abelian lattice gauge theory with tensor networks, Quantum 1, 9 (2017).
  * Bañuls  _et al._ [2017] M. C. Bañuls, K. Cichy, J. I. Cirac, K. Jansen, and S. Kühn, Efficient basis formulation for ($1+1$)-dimensional su(2) lattice gauge theory: Spectral calculations with matrix product states, Phys. Rev. X 7, 041046 (2017).
  * Cataldi  _et al._ [2023] G. Cataldi, G. Magnifico, P. Silvi, and S. Montangero, (2+1)d su(2) yang-mills lattice gauge theory at finite density via tensor networks (2023), arXiv:2307.09396 .
  * Rigobello  _et al._ [2023] M. Rigobello, G. Magnifico, P. Silvi, and S. Montangero, Hadrons in (1+1)d hamiltonian hardcore lattice qcd (2023), arXiv:2308.04488 .
  * Hayata  _et al._ [2023] T. Hayata, Y. Hidaka, and K. Nishimura, Dense $\textrm{QCD}_{2}$ with matrix product states (2023), arXiv:2311.11643 .
  * Zhou  _et al._ [2022] Z.-Y. Zhou, G.-X. Su, J. C. Halimeh, R. Ott, H. Sun, P. Hauke, B. Yang, Z.-S. Yuan, J. Berges, and J.-W. Pan, Thermalization dynamics of a gauge theory on a quantum simulator, Science 377, 311 (2022).
  * Paulson  _et al._ [2021] D. Paulson, L. Dellantonio, J. F. Haase, A. Celi, A. Kan, A. Jena, C. Kokail, R. van Bijnen, K. Jansen, P. Zoller, and C. A. Muschik, Simulating 2d effects in lattice gauge theories on a quantum computer, PRX Quantum 2, 030334 (2021).
  * Angelides  _et al._ [2023b] T. Angelides, P. Naredi, A. Crippa, K. Jansen, S. Kühn, I. Tavernelli, and D. S. Wang, First-order phase transition of the schwinger model with a quantum computer (2023b), arXiv:2312.12831 .
  * Chai  _et al._ [2023] Y. Chai, A. Crippa, K. Jansen, S. Kühn, V. R. Pascuzzi, F. Tacchino, and I. Tavernelli, Entanglement production from scattering of fermionic wave packets: a quantum computing approach (2023), arXiv:2312.02272 [quant-ph] .
  * Davoudi  _et al._ [2024] Z. Davoudi, C.-C. Hsieh, and S. V. Kadam, Scattering wave packets of hadrons in gauge theories: Preparation on a quantum computer (2024), arXiv:2402.00840 [quant-ph] .
  * Ciavarella  _et al._ [2021] A. Ciavarella, N. Klco, and M. J. Savage, Trailhead for quantum simulation of su(3) yang-mills lattice gauge theory in the local multiplet basis, Phys. Rev. D 103, 094501 (2021).
  * Farrell  _et al._ [2023] R. C. Farrell, I. A. Chernyshev, S. J. M. Powell, N. A. Zemlevskiy, M. Illa, and M. J. Savage, Preparations for quantum simulations of quantum chromodynamics in $1+1$ dimensions. ii. single-baryon $\beta$-decay in real time, Phys. Rev. D 107, 054513 (2023).
  * Atas  _et al._ [2023] Y. Y. Atas, J. F. Haase, J. Zhang, V. Wei, S. M.-L. Pfaendler, R. Lewis, and C. A. Muschik, Simulating one-dimensional quantum chromodynamics on a quantum computer: Real-time evolutions of tetra- and pentaquarks, Phys. Rev. Res. 5, 033184 (2023).
  * Ciavarella [2023] A. N. Ciavarella, Quantum simulation of lattice qcd with improved hamiltonians, Phys. Rev. D 108, 094513 (2023).
  * Kruckenhauser  _et al._ [2022] A. Kruckenhauser, R. van Bijnen, T. V. Zache, M. Di Liberto, and P. Zoller, High-dimensional so (4)-symmetric rydberg manifolds for quantum simulation, Quantum Science and Technology 8, 015020 (2022).
  * Cohen and Thompson [2021] S. R. Cohen and J. D. Thompson, Quantum computing with circular rydberg atoms, PRX Quantum 2, 030322 (2021).
  * Chi  _et al._ [2022] Y. Chi, J. Huang, Z. Zhang, J. Mao, Z. Zhou, X. Chen, C. Zhai, J. Bao, T. Dai, H. Yuan, _et al._ , A programmable qudit-based quantum processor, Nature communications 13, 1166 (2022).
  * Kasper  _et al._ [2021] V. Kasper, D. González-Cuadra, A. Hegde, A. Xia, A. Dauphin, F. Huber, E. Tiemann, M. Lewenstein, F. Jendrzejewski, and P. Hauke, Universal quantum computation and quantum error correction with ultracold atomic mixtures, Quantum Science and Technology 7, 015008 (2021).
  * Low  _et al._ [2023] P. J. Low, B. White, and C. Senko, Control and readout of a 13-level trapped ion qudit, (2023), arXiv:2306.03340 .
  * González-Cuadra  _et al._ [2022] D. González-Cuadra, T. V. Zache, J. Carrasco, B. Kraus, and P. Zoller, Hardware efficient quantum simulation of non-abelian gauge theories with qudits on rydberg platforms, Phys. Rev. Lett. 129, 160501 (2022).
  * Zache  _et al._ [2023] T. V. Zache, D. González-Cuadra, and P. Zoller, Fermion-qudit quantum processors for simulating lattice gauge theories with matter, (2023), arXiv:2303.08683 .
  * Popov  _et al._ [2023] P. P. Popov, M. Meth, M. Lewenstein, P. Hauke, M. Ringbauer, E. Zohar, and V. Kasper, Variational quantum simulation of u (1) lattice gauge theories with qudit systems, (2023), arXiv:2307.15173 .
  * Meth  _et al._ [2023] M. Meth, J. F. Haase, J. Zhang, C. Edmunds, L. Postler, A. Steiner, A. J. Jena, L. Dellantonio, R. Blatt, P. Zoller, T. Monz, P. Schindler, C. Muschik, and M. Ringbauer, Simulating 2D lattice gauge theories on a qudit quantum computer, (2023), arXiv:2310.12110 .
  * Zalivako  _et al._ [2024] I. V. Zalivako, A. S. Nikolaeva, A. S. Borisenko, A. E. Korolkov, P. L. Sidorov, K. P. Galstyan, N. V. Semenin, V. N. Smirnov, M. A. Aksenov, K. M. Makushin, _et al._ , Towards multiqudit quantum processor based on a ${}^{171}\textrm{Yb}^{+}$ ion string: Realizing basic quantum algorithms, (2024), arXiv:2402.03121 .
  * Low  _et al._ [2020] P. J. Low, B. M. White, A. A. Cox, M. L. Day, and C. Senko, Practical trapped-ion protocols for universal qudit-based quantum computing, Physical Review Research 2, 033128 (2020).
  * Susskind [1977] L. Susskind, Lattice fermions, Phys. Rev. D 16, 3031 (1977).
  * Zohar and Burrello [2015] E. Zohar and M. Burrello, Formulation of lattice gauge theories for quantum simulations, Phys. Rev. D 91, 054506 (2015).
  * Chandrasekharan and Wiese [1997] S. Chandrasekharan and U.-J. Wiese, Quantum link models: A discrete approach to gauge theories, Nuclear Physics B 492, 455 (1997).
  * Silvi  _et al._ [2019b] P. Silvi, Y. Sauer, F. Tschirsich, and S. Montangero, Tensor network simulation of an su(3) lattice gauge theory in 1d, Phys. Rev. D 100, 074512 (2019b).
  * Zohar and Cirac [2019] E. Zohar and J. I. Cirac, Removing staggered fermionic matter in $u(n)$ and $su(n)$ lattice gauge theories, Phys. Rev. D 99, 114511 (2019).
  * Mølmer and Sørensen [1999] K. Mølmer and A. Sørensen, Multiparticle entanglement of hot trapped ions, Phys. Rev. Lett. 82, 1835 (1999).
  * Roos  _et al._ [1999] C. Roos, T. Zeiger, H. Rohde, H. C. Nägerl, J. Eschner, D. Leibfried, F. Schmidt-Kaler, and R. Blatt, Quantum state engineering on an optical transition and decoherence in a paul trap, Phys. Rev. Lett. 83, 4713 (1999).
  * Schindler  _et al._ [2013] P. Schindler, D. Nigg, T. Monz, J. T. Barreiro, E. Martinez, S. X. Wang, S. Quint, M. F. Brandl, V. Nebendahl, C. F. Roos, M. Chwalla, M. Hennrich, and R. Blatt, A quantum information processor with trapped ions, New Journal of Physics 15, 123012 (2013).
  * Häffner  _et al._ [2003] H. Häffner, S. Gulde, M. Riebe, G. Lancaster, C. Becher, J. Eschner, F. Schmidt-Kaler, and R. Blatt, Precision measurement and compensation of optical stark shifts for an ion-trap quantum processor, Phys. Rev. Lett. 90, 143602 (2003).
  * Kirchmair  _et al._ [2009] G. Kirchmair, J. Benhelm, F. Zähringer, R. Gerritsma, C. F. Roos, and R. Blatt, Deterministic entanglement of ions in thermal states of motion, New Journal of Physics 11, 023002 (2009).
  * Ruster  _et al._ [2016] T. Ruster, C. T. Schmiegelow, H. Kaufmann, C. Warschburger, F. Schmidt-Kaler, and U. G. Poschinger, A long-lived zeeman trapped-ion qubit, Applied Physics B 122, 254 (2016).
  * Moses [2023] S. A. e. a. Moses, A race-track trapped-ion quantum processor, Phys. Rev. X 13, 041052 (2023).
  * Gerster  _et al._ [2022] L. Gerster, F. Martínez-García, P. Hrmo, M. W. van Mourik, B. Wilhelm, D. Vodola, M. Müller, R. Blatt, P. Schindler, and T. Monz, Experimental bayesian calibration of trapped-ion entangling operations, PRX Quantum 3, 020350 (2022).
  * Klco  _et al._ [2020] N. Klco, M. J. Savage, and J. R. Stryker, Su(2) non-abelian gauge field theory in one dimension on digital quantum computers, Phys. Rev. D 101, 074512 (2020).
  * Liu  _et al._ [2023] H. Liu, T. Bhattacharya, S. Chandrasekharan, and R. Gupta, Phases of 2d massless qcd with qubit regularization, (2023), arXiv:2312.17734 .
  * Bauer  _et al._ [2023] C. W. Bauer, I. D’Andrea, M. Freytsis, and D. M. Grabowska, A new basis for hamiltonian su (2) simulations, (2023), arXiv:2307.11829 .
  * Mathis  _et al._ [2020] S. V. Mathis, G. Mazzola, and I. Tavernelli, Toward scalable simulations of lattice gauge theories on quantum computers, Phys. Rev. D 102, 094501 (2020).
  * Fishman  _et al._ [2022a] M. Fishman, S. R. White, and E. M. Stoudenmire, The ITensor Software Library for Tensor Network Calculations, SciPost Phys. Codebases , 4 (2022a).
  * Fishman  _et al._ [2022b] M. Fishman, S. R. White, and E. M. Stoudenmire, Codebase release 0.3 for ITensor, SciPost Phys. Codebases , 4 (2022b).