#  Quantum computational resources for lattice QCD in the strong-coupling
limit

Michael Fromm<EMAIL_ADDRESS>Institut für Theoretische Physik,
Goethe-Universität Frankfurt
Max-von-Laue-Str. 1, 60438 Frankfurt am Main, Germany  Lucas Katschke
<EMAIL_ADDRESS>Institut für Theoretische Physik, Goethe-
Universität Frankfurt
Max-von-Laue-Str. 1, 60438 Frankfurt am Main, Germany  Owe Philipsen
<EMAIL_ADDRESS>Institut für Theoretische Physik, Goethe-
Universität Frankfurt
Max-von-Laue-Str. 1, 60438 Frankfurt am Main, Germany  John von Neumann
Institute for Computing (NIC) at GSI
Planckstr. 1, 64291 Darmstadt, Germany  Wolfgang Unger<EMAIL_ADDRESS>bielefeld.de Fakultät für Physik, Bielefeld University
33615 Bielefeld, Germany

###### Abstract

We consider the strong coupling limit of lattice QCD with massless staggered
quarks and study the resource requirements for quantum simulating the theory
in its Hamiltonian formulation. The bosonic Hilbert space of the color-singlet
degrees of freedom grows quickly with the number of quark flavors, making it a
suitable testing ground for resource considerations across different
platforms. In particular, in addition to the standard model of computation
with qubits, we consider mapping the theory to qudits $(d>2)$ and qumodes, as
used on trapped-ion systems and photonic devices, respectively.

Hamiltonian lattice gauge theory, Quantum computing

##  I Introduction

It was a key observation that universal quantum computers offer the potential
to simulate aspects of quantum systems efficiently that are intractable by
classical machines Lloyd (1996), provided the system’s interaction remains
local. As many of these are relevant for High Energy Physics, it is not
surprising that quantum computation quickly became the focus of a substantial
theoretical effort in the field, where we point to the whitepapers Bauer  _et
al._ (2022); Di Meglio  _et al._ (2023) that summarize the progress and lay
out the roadmap. In this context, lattice gauge theory Wilson (1974) (LGT)
furnishes a non-perturbative approach to many problems at hand. Even though
the Hamiltonian formulation was derived very early Kogut and Susskind (1975),
its application to universal quantum computation was realized and examined
under this aspect much later Byrnes and Yamamoto (2005). By simple
considerations it became clear that contrary to the use of LGT on classical
computers, which is often hampered by the fermionic sign problem, the mapping
of the bosonic degrees of freedom (dof) to the quantum computational dof
requires truncation. However, the approximation involved is _controllable_ and
systematic error bounds e.g. for Hamiltonian time evolution can be derived
Tong  _et al._ (2021). Given this solid theoretical ground, the development of
efficient formulations for non-Abelian LGT is hence of great interest Davoudi
_et al._ (2020) and ranges from discrete subgroups Assi and Lamm (2024) for
pure gauge theory to formulations in terms of gauge invariant dof s.a. the
loop-string-hadron formulation Davoudi  _et al._ (2023) for theories with
matter.

In this work we investigate another aspect of the mapping of lattice gauge
theory and estimate the resource requirements if we vary the standard
computational unit. In particular, besides two-level qubits we consider multi-
level qudits with a state space of $d>2$ Wang  _et al._ (2020) and quantum
modes (qumodes) Weedbrook  _et al._ (2012) as used for continuous-variable
quantum computation Lloyd and Braunstein (1999) on photonic devices.

Testing ground for our study will be the strong-coupling limit of lattice QCD
as re-examined recently in Fromm  _et al._ (2023). Actual interest in this
theory arose with the work of Rossi and Wolff (1984); Karsch and Mütter
(1989), back in the early days of Euclidean LGT, mainly motivated by the fact
that the fermionic sign problem of the theory remained tractable to some
degree. As a consequence the phase diagram of the theory with $N_{f}=1$
staggered quark flavor could be obtained de Forcrand and Fromm (2010), soon
followed by the same study including gauge corrections de Forcrand  _et al._
(2014).

By taking the continuous-time limit of the Euclidean LGT, the Hamiltonian for
the theory with varying number of quark flavors could be constructed Klegrewe
and Unger (2020); Unger and Pattanaik (2022a, b). It is this last development
which makes the theory interesting from the perspective of quantum simulation:
The resulting Hilbert space is a discrete state-space of physical, color-
singlet dof (mesons and baryons) whose interaction is controlled by a local,
sparse nearest-neighbor Hamiltonian. However, the Hilbert space dimension of
the theory grows quickly with increasing number of staggered quark flavors
$N_{f}$ and is thus quantum computationally well beyond the hardware limits of
the NISQ era. For this reason we study the theory and its resource
requirements after a suitable map onto different computational units (qubit,
qudit or qumode). Before we lay out the results of this analysis in Sect.III,
we present a short account of the derivation of the strong coupling theory in
the Hamiltonian formulation in the following Sect.II.

##  II Background

In Ref.Klegrewe and Unger (2020) the dual formulation of strong coupling
Lattice QCD with $N_{f}=1$ flavors of staggered quarks was derived in the
limit of continuous Euclidean time by letting the lattice extent in time
direction $N_{\tau}\rightarrow\infty$, while keeping the temperature in
spatial lattice units $aT\sim\gamma^{2}/N_{t}$ fixed, where $\gamma$
represents the bare anisotropy parameter. The resulting partition function at
temperature $T$ and baryon chemical potential $\mu_{B}$,

$\displaystyle\mathcal{Z}_{\rm CT}(aT,a\mu_{B})$
$\displaystyle=\operatorname{Tr}_{\mathbbm{H}_{\mathfrak{h}}}\left[e^{(-\hat{\mathcal{H}}+\hat{\mathcal{N}}a\mu_{B})/aT}\right],$
(1)

gives then rise to a Hamiltonian and number operator

$\displaystyle\hat{\mathcal{H}}$ $\displaystyle=-\frac{1}{2}\sum_{\langle
x,y\rangle}\left(\hat{J}^{+}_{x}\hat{J}^{-}_{y}+\hat{J}^{-}_{x}\hat{J}^{+}_{y}\right),$
$\displaystyle\hat{\mathcal{N}}$ $\displaystyle=\sum_{x}\hat{\omega}_{x}.$ (2)

with matrix representations

$\displaystyle\hat{J}^{+}_{x}$
$\displaystyle=\left(\begin{array}[]{cccc|cc}0&0&0&0&&\\\ v_{L}&0&0&0&&\\\
0&v_{T}&0&0&&\\\ 0&0&v_{L}&0&&\\\ \hline\cr&&&&0&0\\\ &&&&0&0\\\
\end{array}\right),\,\hat{\omega}=\left(\begin{array}[]{cccc|cc}0&0&0&0&&\\\
0&0&0&0&&\\\ 0&0&0&0&&\\\ 0&0&0&0&&\\\ \hline\cr&&&&1&0\\\
&&&&0&-1\end{array}\right).$ (15)

The trace in Eq. (1) is over the Hilbert space
${\mathbbm{H}_{\mathfrak{h}}}=\bigotimes_{x\in\Sigma}{\mathbbm{H}_{\mathfrak{h}}}_{x}$
where ${\mathbbm{H}_{\mathfrak{h}}}_{x}$ is six dimensional, with states
$|\mathfrak{h}\rangle=|\mathfrak{m},\mathfrak{b}\rangle\in{\mathbbm{H}_{\mathfrak{h}}}_{x}$
and a basis $|\mathfrak{h}_{i}\rangle$ of six hadronic states (four mesonic
and two baronic) per spatial lattice site,

$\displaystyle|\mathfrak{h}_{i}\rangle\in\\{0,\pi,2\pi,3\pi,B^{+},B^{-}\\}\,.$
(16)

With the mesonic occupation number $\mathfrak{m}=0,\pi,2\pi,3\pi$ being also
bounded from above, the theory has a particle-hole symmetry. If we let
$\mathfrak{s}=\mathfrak{m}-3/2$, we see that the mesonic part of the
Hamiltonian in Eq.(2) actually corresponds to the quantum Heisenberg model in
spin $j=N_{c}/2=3/2$ representation with states $|j,\mathfrak{s}\rangle$.111We
note that the antiferromagnetic Heisenberg Hamiltonian also emerges as an
effective theory of Lattice QCD in the strong coupling limit in lowest order
of Hamiltonian perturbation theory Smit (1980).

Also, owing to the limit of infinite bare gauge coupling,
$\beta=2N_{c}/g^{2}=0$, baryons become static with trivial dynamics in the
continuous time limit - a fact, that allows for importance sampling of the
theory at finite baryon density _without_ a sign problem. As a consequence,
$\hat{\mathcal{H}}$ and $\hat{\mathcal{N}}$ in Eq.(2) have a block-diagonal
structure and the states $\ket{\mathfrak{h}}_{x}$ of the local Hilbert space
are a direct sum of mesonic and baryonic occupational states,
$\ket{\mathfrak{h}}_{x}=\ket{\mathfrak{m}}_{x}\oplus\ket{\mathfrak{b}}_{x}$.

In Unger and Pattanaik (2022a, b); Unger (2024) the theory was extended to
$N_{f}=2$ flavors of staggered quarks. The partition function is now given by

$\displaystyle Z_{\rm CT}(aT,a\mu_{B},a\mu_{I})$
$\displaystyle=\operatorname{Tr}_{\mathbbm{H}_{\mathfrak{h}}}\left[e^{(-\hat{\mathcal{H}}+\hat{\mathcal{N}}_{B}a\mu_{B}+\hat{\mathcal{N}}_{I}a\mu_{I})/aT}\right]$
(17)

where the Hamiltonian is a sum of four contributions, one for each “pion”
$\pi\in\\{\pi^{+},\pi^{-},\pi_{u},\pi_{d}\\}$:

$\displaystyle\hat{\mathcal{H}}$ $\displaystyle=-\frac{1}{2}\sum_{\langle
x,y\rangle}\sum_{Q_{i}\in\\{\begin{subarray}{c}\pi^{+},\pi^{-},\\\
\pi_{U},\pi_{D}\end{subarray}\\}}\left({\hat{J}_{Q_{i},x}}^{+}{\hat{J}_{Q_{i},y}}^{-}+{\hat{J}_{Q_{i},x}}^{-}{\hat{J}_{Q_{i},y}}^{+}\right).$
(18)

In addition to the baryon number operator $\hat{\mathcal{N}}_{B}={\rm
diag}(-2,-1,\ldots 1,2)$, the isospin number operator
$\hat{\mathcal{N}}_{I}={\rm
diag}\left(0,-\frac{3}{2},\ldots\frac{3}{2},0\right)$ has been introduced and
coupled to the isospin chemical potential $a\mu_{I}$. The local Hilbert space
${\mathbbm{H}_{\mathfrak{h}}}_{x}$ now has a basis of 92 distinct hadronic
states $|\mathfrak{h}\rangle$. The particle-hole symmetry of the one-flavor
case with quantum number $\mathfrak{s}$ generalizes to
$\mathfrak{s}=\mathfrak{m}-\frac{N_{c}}{2}(N_{f}-|B|)$, in addition the states
can be labelled by the local baryon number $B$ and isospin $I$ as shown in
Table 1. Once again, as the baryon number is conserved and baryons remain
static, the ladder operators $\hat{J}_{Q_{i},x}^{+/-}$ in Eq.(18) have block-
diagonal structure and decompose into sectors of fixed local baryon number
$B=-2,\ldots,2$. In each sector the ladder operators generate a corresponding
state multiplet which is visible in Fig.1-Fig.1. The action of the ladder
operators $J^{+}_{Q_{i}}$ on a state $|\mathfrak{h}\rangle$ is indicated by
yellow, blue, green and red arrows for the conserved currents
$Q_{i}=\pi_{u},\pi_{d},\pi^{+}$ and $\pi^{-}$, respectively. Apart from the
state number, a state is labelled by its mesonic and baryonic occupation
number where we point to Unger and Pattanaik (2022a, b) for a full account of
the derivation of the matrix elements and the state definition. For our
purpose, it suffices to know that the $\hat{J}^{+}_{Q_{i}}$ once again
generate a representation of a $SU(2)$ algebra, in particular

$\displaystyle\frac{N_{c}}{2}[\hat{J}^{+}_{Q_{i}},\hat{J}^{-}_{Q_{j}}]=\hat{J}^{3}_{Q_{i}}\delta_{ij}\,.$
(19)

This time however the nature of the representation is more intricate. With the
particular labelling chosen in Fig.1-Fig.1, where isospin $I$ increases
horizontally from left to right and $\mathfrak{s}$ increases vertically, c.f.
Table 1, we see that the representation of the $J^{+}_{\pi_{u}}$ and
$J^{+}_{\pi_{d}}$, which are isospin conserving, is again reducible and hence
becomes block-diagonal in each isospin sector, where it in turn forms a
product representation.

Figure 1: Hadronic states in the sectors with baryon number $B=0$ $(a)$, $B=1,2$ $(b)$ and $B=-1,-2$ $(c)$. The action of the ladder operators $J^{+}_{Q_{i}}$ on a state is indicated by yellow, blue, green and red arrows for the conserved currents $Q_{i}=\pi_{u},\pi_{d},\pi^{+}$ and $\pi^{-}$, respectively. In each sector, Isospin $I$ increases horizontally from left to right, the mesonic occupation number $\mathfrak{s}$ increases vertically, bottom to top. $B$ | $I$ | $\mathfrak{s}=\mathfrak{m}-\frac{3}{2}(2-|B|)$ | $\Sigma$
---|---|---|---
|  | $-3$ | $-\frac{5}{2}$ | $-2$ | $-\frac{3}{2}$ | $-1$ | $-\frac{1}{2}$ | $\,0\,$ | $+\frac{1}{2}$ | $+1$ | $+\frac{3}{2}$ | $+2$ | $+\frac{5}{2}$ | $+3$ |
-2 | 0 |  |  |  |  |  |  | 1 |  |  |  |  |  |  | 1
-1 | $-\frac{3}{2}$ |  |  |  | 1 |  | 1 |  | 1 |  | 1 |  |  |  | 4
-1 | $-\frac{1}{2}$ |  |  |  | 1 |  | 2 |  | 2 |  | 1 |  |  |  | 6
-1 | $+\frac{1}{2}$ |  |  |  | 1 |  | 2 |  | 2 |  | 1 |  |  |  | 6
-1 | $+\frac{3}{2}$ |  |  |  | 1 |  | 1 |  | 1 |  | 1 |  |  |  | 4
0 | $-3$ |  |  |  |  |  |  | 1 |  |  |  |  |  |  | 1
0 | $-2$ |  |  |  |  | 1 |  | 2 |  | 1 |  |  |  |  | 4
0 | $-1$ |  |  | 1 |  | 2 |  | 4 |  | 2 |  | 1 |  |  | 10
0 | $0$ | 1 |  | 2 |  | 4 |  | 6 |  | 4 |  | 2 |  | 1 | 20
0 | $+1$ |  |  | 1 |  | 2 |  | 4 |  | 2 |  | 1 |  |  | 10
0 | $+2$ |  |  |  |  | 1 |  | 2 |  | 1 |  |  |  |  | 4
0 | $+3$ |  |  |  |  |  |  | 1 |  |  |  |  |  |  | 1
1 | $-\frac{3}{2}$ |  |  |  | 1 |  | 1 |  | 1 |  | 1 |  |  |  | 4
1 | $-\frac{1}{2}$ |  |  |  | 1 |  | 2 |  | 2 |  | 1 |  |  |  | 6
1 | $+\frac{1}{2}$ |  |  |  | 1 |  | 2 |  | 2 |  | 1 |  |  |  | 6
1 | $+\frac{3}{2}$ |  |  |  | 1 |  | 1 |  | 1 |  | 1 |  |  |  | 4
2 | 0 |  |  |  |  |  |  | 1 |  |  |  |  |  |  | 1
$\Sigma$ |  | 1 | 0 | 4 | 8 | 10 | 12 | 22 | 12 | 10 | 8 | 4 | 0 | 1 | 92
Table 1:  State multiplicities of the local state space for the two-flavor
theory in the Hamiltonian formulation, totalling $92$ states at site $x$. The
states are labelled by baryon number $B$, isospin number $I$ and symmetrized
meson occupation number
$\mathfrak{s}=\mathfrak{m}-\frac{N_{\text{c}}}{2}(N_{\text{f}}-|B|)$.

##  III Results

###  III.1 $N_{f}=1$

####  III.1.1 General Remarks

Figure 2:  Schematic of a single Trotter step of step size $\delta t$,
connecting the $nn$-sites $x$ and $y$ with local states
$|\mathfrak{h}\rangle_{x}$ and $|\mathfrak{h}\rangle_{y}$.

The operators appearing in Eq.(15) have block-diagonal structure: The meson
and baryon sector do not mix when acting with the Hamiltonian on the product
space
${\mathbbm{H}_{\mathfrak{h}}}_{x}\otimes{\mathbbm{H}_{\mathfrak{h}}}_{y}$. If
we extend the $|B|>0$ sector also to 4 dimensions, we can write the one-flavor
Hamiltonian in terms of $|B|$-projectors,

$\displaystyle\hat{\mathcal{H}}=-\frac{1}{2}\sum_{\langle\vec{x},\vec{y}\rangle}|0\rangle\langle
0|_{x}|0\rangle\langle
0|_{y}\left(\hat{J}^{+}_{\vec{x}}\hat{J}^{-}_{\vec{y}}+\hat{J}^{-}_{\vec{x}}\hat{J}^{+}_{\vec{y}}\right)$
$\displaystyle+a\mu_{B}\sum_{\vec{x}}|1\rangle\langle 1|_{x}\hat{\omega}_{x}.$
(20)

We used that property to derive the gate set for $N_{f}=1$ in terms of
controlled gates where locally this amounts to an ancillary bit, controlling
the baryonic sector Fromm  _et al._ (2023).

####  III.1.2 Mapping to Qudits

Using a $d$-level systems with $d>2$ (qudit) as quantum computational unit,
the dimension of the local Hilbert space in Eq.(16) suggests the use of qudits
with $d=6$. Rewriting the Hamiltonian in terms of
$\hat{J}^{1}=\frac{\sqrt{3}}{2}(J^{+}+J^{-})$ and
$\hat{J}^{2}=\frac{\sqrt{3}}{2i}(\hat{J}^{+}-\hat{J}^{-})$, the $\hat{J}^{i}$
can be decomposed into generators $\mathcal{X}_{ij},\mathcal{Y}_{ij}$ of
unitary Givens rotations $G_{jk}(\theta,\phi)$, with

$\displaystyle
G_{jk}(\theta,\phi)=e^{-i\theta(\cos{\phi}\mathcal{X}_{jk}-\sin{\phi}\mathcal{Y}_{jk})}\,,$
(21)

where O’Leary  _et al._ (2006)

$\displaystyle\mathcal{X}_{kl}=|k\rangle\langle l|+|l\rangle\langle
k|,\,\mathcal{Y}_{kl}=-i(|k\rangle\langle l|-|l\rangle\langle k|)\,.$ (22)

One obtains

$\displaystyle J^{1}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{3}}{2}(\mathcal{X}_{01}+\frac{2}{\sqrt{3}}\mathcal{X}_{12}+\mathcal{X}_{23})$
$\displaystyle J^{2}$ $\displaystyle=$
$\displaystyle-\frac{\sqrt{3}}{2}(\mathcal{Y}_{01}+\frac{2}{\sqrt{3}}\mathcal{Y}_{12}+\mathcal{Y}_{23})\,.$
(23)

If we trotterize the time evolution for a pair of neighboring mesonic sites
$x,y$ with the hopping term in Eq.(2) according to

$\displaystyle e^{-i\delta
tH_{x,y}}=\prod_{i}U_{\mathcal{XX}_{i}}U_{\mathcal{YY}_{i}}+\mathcal{O}(\delta
t^{2})\,$ (24)

we obtain a decomposition into 9 two-qudit entangling operations of the
generic form
$U_{\mathcal{XX}}=e^{-i\alpha\mathcal{X}^{x}_{kl}\otimes\mathcal{X}^{y}_{mn}}$
and
$U_{\mathcal{YY}}=e^{-i\beta\mathcal{Y}^{x}_{kl}\otimes\mathcal{Y}^{y}_{mn}}$,
respectively. These elementary operations can be realized via phase-
compensated Mølmer Sørensen (MS) gates for trapped-ion qudits Molmer and
Sorensen (1998); Ringbauer  _et al._ (2021). In Low  _et al._ (2019) a
generalization of the MS gate was proposed that allows for the simultaneous
driving of multiple transitions between qudit levels of the involved
neighboring qudits in an entangling operation. In other words, a linear
combination of two-level transitions s.a. Eqs.(23) can be realized by just one
operation. As a result, the hopping term can be implemented via two two-qudit
entangling-operations

$\displaystyle e^{-i\delta tH_{x,y}}=e^{-i\delta tJ^{1}_{x}\otimes
J^{1}_{y}}e^{-i\delta tJ^{2}_{x}\otimes J^{2}_{y}}+\mathcal{O}(\delta
t^{2}),\,$ (25)

thus reducing the gate depth to 2 for the $nn$-term. This substantial
reduction in entangling operations comes at the cost of increased complexity
in experimental control, where we refer to the analysis carried out in Calajò
_et al._ (2024) to address the experimental challenges.222See also Illa  _et
al._ (2024) for a recent resource analysis of digital quantum simulations of
(1+1)D Lattice QCD on $d=8$ trapped ion qudits.

In contrast, the number operator $\hat{\mathcal{N}}$ in Eq.(2) just
corresponds to a single, diagonal qudit operation. Picking as a basis for
diagonal gates the $\mathcal{Z}_{jk}$ with

$\displaystyle\mathcal{Z}_{jl}=|j\rangle\langle j|-|l\rangle\langle l|$ (26)

that, together with $\mathcal{X}_{jl}$ and $\mathcal{Y}_{jl}$, generate a
$su(2)$ subalgebra of $SU(d=6)$, the local term in Eq.(2) then generates an
one-qudit operation of the type $e^{-it\mu_{B}\mathcal{Z}_{45}}$ for each
qudit during time evolution. As this corresponds to a single diagonal gate
without Trotterization, the main cost of the time evolution will be given by
Eq.(25).

####  III.1.3 Mapping to Qumodes

As a third alternative, we can encode our system in the degrees of freedom of
quantum harmonic oscillators, i.e. qumodes. The basic operator set then
consists of the quadrature variables $(\hat{x_{i}},\hat{p_{i}})$ or,
alternatively, Fock-space creation and annihilation operators,
$\hat{a}_{i},\hat{a}_{i}^{\dagger}$, with the usual relations

$\displaystyle\hat{x}_{i}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(\hat{a}_{i}+\hat{a}_{i}^{\dagger}),\,\hat{p}_{i}=\frac{1}{2i}(\hat{a}_{i}-\hat{a}_{i}^{\dagger})$
(27) $\displaystyle\hat{a}_{i}$ $\displaystyle=$ $\displaystyle
x_{i}+ip_{i},\,a_{i}^{\dagger}=x_{i}-ip_{i}$ (28)

and

$\displaystyle[\hat{a}_{i},\hat{a}_{j}]=[\hat{a}_{i}^{\dagger},\hat{a}_{j}^{\dagger}]=0,[\hat{a}_{i},\hat{a}^{\dagger}_{j}]=\delta_{ij},{[\hat{x}_{i},\hat{p}_{j}]}=\frac{i}{2}\delta_{ij}\,.$

The $j=\frac{N_{c}}{2}$ representation of the $J^{+}$ and $J^{-}$ (non-trivial
block in Eq.(15)) can be obtained via the Jordan-Schwinger-map Schwinger
(1952); Annabestani  _et al._ (2020)

$\displaystyle\hat{J}_{x}^{+}=\hat{a}_{1,x}^{\dagger}\hat{a}_{2,x},\,\hat{J}_{x}^{-}=\hat{a}_{2,x}^{\dagger}\hat{a}_{1,x},$
(29)

s.t.

$\displaystyle|j,m\rangle=\frac{(\hat{a}_{1}^{\dagger})^{j+m}(\hat{a}_{2}^{\dagger})^{j-m}}{\sqrt{(j+m)!(j-m)!}}|0\rangle\,.$
(30)

The local number operators
$\hat{n}_{i}=\hat{a}_{i}^{\dagger}\hat{a}_{i},\,i=1,2$ then have eigenvalues
$n_{1/2}$ which are constrained to $n_{1}+n_{2}=2j=N_{c}$. The mesonic hopping
term of the Hamiltonian Eq.(2) hence reads

$\displaystyle\mathcal{H}_{meson}=-\frac{1}{2}\sum_{\langle{x},{y}\rangle}(\hat{a}^{\dagger}_{1,x}\hat{a}_{2,x}\hat{a}^{\dagger}_{2,y}\hat{a}_{1,y})+h.c.\,.$
(31)

Expressed in the quadrature field operators Eq.(27) this expression becomes

$\displaystyle\mathcal{H}_{meson}=-\frac{1}{2}\sum_{\langle{x},{y}\rangle}\sum_{\\{i_{k}\\}}(\hat{q}^{i_{1}}_{1,x}\hat{q}^{i_{2}}_{2,x}\hat{q}^{i_{3}}_{2,y}\hat{q}^{i_{4}}_{1,y})s_{\\{i_{1},\ldots,i_{4}\\}},$
(32)

with the shorthand $(\hat{q}^{1},\hat{q}^{2})=(\hat{x},\hat{p})$ and signs
$s_{\\{i_{1},\ldots,i_{4}\\}}$333$s_{\\{2,1,1,2\\}}=s_{\\{1,2,2,1\\}}=-1$,
otherwise $s_{\\{i_{1},\ldots,i_{4}\\}}=1$.. In the following we assume to
have access to the qumode instruction set

$\displaystyle\\{e^{is\hat{x}_{k}},e^{is\hat{a}_{k}^{\dagger}\hat{a}_{k}},e^{is\hat{x}_{k}^{2}},e^{is\hat{x}_{k}^{3}},e^{is\hat{x}_{k}\hat{x}_{j}}\\}$
$\displaystyle\equiv\\{Z(s),R(s),P(s),V(s),CZ(s)\\}$ (33)

In particular, the Fourier gate $F\equiv
e^{i\frac{\pi}{2}\hat{a}_{k}^{\dagger}\hat{a}_{k}}$ can be used for unitary
conjugation of powers of quadrature operators, e.g.
$g(\hat{p}_{k})=F_{k}g(\hat{x}_{k})F_{k}^{\dagger}$ for polynomials $g$. This
allows us to rewrite each of the eight terms in Eq.(32) in the form
$\hat{x}_{1,x}\hat{x}_{2,x}\hat{x}_{1,y}\hat{x}_{2,y}$. Through exact gate
decompositions (Appendix A.1) trotterized time evolution by a single term
$e^{i\delta t\hat{x}_{1,x}\hat{x}_{2,x}\hat{x}_{1,y}\hat{x}_{2,y}}$ is found
to be costly, involving approximately $8\times 65$ elementary operations of
the instruction set Eq.(33), which reduces to $8\times 11$ elementary
operations, depending on the availability of a quartic gate
$Q(s)=e^{is\hat{x}_{k}^{4}}$. In total, a first-order Trotter step of a single
mesonic $nn$-pair, Eq.(50), results in a gate count of $64Q,147CZ,149F$ and a
corresponding gate depth of 348 (see Appendix A.1), now under the assumption
that the quartic gate is part of the gate set.

Similarily, the baryonic term of Eq.(2) can be mapped to

$\displaystyle\mathcal{H}_{baryon}=\sum_{x}(\hat{n}_{1,x}-\hat{n}_{2,x})=\sum_{x,i}(\hat{x}^{2}_{i,x}+\hat{p}_{i,x}^{2})$
(34)

now applied to Fock-states with the constraint $n_{1}+n_{2}=1$. The
exponentiated term can be implemented using two rotation gates
$R(\mu_{B}\delta t)=e^{i\mu_{B}\delta t\hat{a}_{k}^{\dagger}\hat{a}_{k}}$ per
baryonic site $x$.

In order to combine the mesonic and baryonic time evolution, we can follow the
strategy of Sect.III.1.1 and control the local sector (mesonic or baryonic)
via an ancillary qubit, yielding a Hamiltonian

$\displaystyle\mathcal{H}=-$ $\displaystyle\frac{1}{2}$
$\displaystyle\sum_{\langle{x},{y}\rangle}|0\rangle\langle
0|_{x}|0\rangle\langle
0|_{y}\sum_{\\{i_{k}\\}}(\hat{q}^{i_{1}}_{1,x}\hat{q}^{i_{2}}_{2,x}\hat{q}^{i_{3}}_{2,y}\hat{q}^{i_{4}}_{1,y})s_{\\{i_{1}\ldots
i_{4}\\}}$ (35) $\displaystyle+$ $\displaystyle
a\mu_{B}\sum_{x}|1\rangle\langle 1|_{x}(\hat{n}_{1,x}-\hat{n}_{2,x})$

Such control implies a mixed qubit-qumode hardware architecture as used in
circuit QED (cQED) Blais  _et al._ (2020).

Alternatively, we could use the fact that either part of the Hamiltonian,
$\mathcal{H}_{meson}$ and $\mathcal{H}_{baryon}$, commutes with the local
number operator $\hat{n}=\hat{n}_{1}+\hat{n}_{2}$ and hence does not change
the local occupation number $n=n_{1}+n_{2}=2j$ of a single Fock-state.
Consequently, we can use the number operator to control the application of the
time evolution operator in the baryonic and mesonic sector _without_ ancillary
qubits, by letting e.g.

$\displaystyle\mathcal{H}_{meson}$
$\displaystyle=-\frac{1}{2}\sum_{\langle{x},{y}\rangle}H_{x,y}\frac{1}{2}(\hat{n}_{x}-1)\otimes\frac{1}{2}(\hat{n}_{y}-1)$
(36) $\displaystyle\mathcal{H}_{baryon}$
$\displaystyle=a\mu_{B}\sum_{x}\mathcal{N}_{x}\frac{1}{2}(3-\hat{n}_{x})\,,$

in a manner reminiscent of the definition of controlled qumode operations s.a.
the controlled Kerr gate Killoran  _et al._ (2019). In the above Eq.(36) we
have implicitly defined the local hopping term $H_{xy}$ and on-site terms
$\mathcal{N}_{x}$ from Eqs.(31) and (34), respectively. We note that the terms
in Eq.(36) then actually correspond to a higher order interaction. In view of
the already considerable cost to trotterize a time evolution with Eq.(31), we
have not pursued this further.

###  III.2 $N_{f}=2$

####  III.2.1 General remarks, mapping to qubits

When mapping the two-flavor theory to qubits, we can use the block-diagonal
shape of the ladder operators, c.f. Fig.1, and rewrite the Hamiltonian Eq.(18)
in terms of projectors of baryonic sectors, in complete analogy to the one-
flavor case, Eq.(20):

$\displaystyle\hat{\mathcal{H}}$ $\displaystyle=-\frac{1}{2}\sum_{\langle
x,y\rangle}\sum_{i,j}|i\rangle\langle i|_{x}|j\rangle\langle
j|_{y}\sum_{Q_{k}}\left({\hat{J}_{Q_{k},x}}^{(i)+}{\hat{J}_{Q_{k},y}}^{(j)-}+{\hat{J}_{Q_{k},x}}^{(i),-}{\hat{J}_{Q_{k},y}}^{(j),+}\right)$
(37)

In the above expression we have expanded each sector’s dimension by "zero
padding" to match the dimensionality of the largest sector ($B=0$) which is
50-dimensional. Formally, the expression Eq.(37) can be derived by noting that
locally for $x\in\Sigma$ we have $J^{+}_{Q_{k},x}=\sum_{i}|i\rangle\langle
i|_{x}\otimes J^{(i)+}_{Q_{k},x}$ where we wrote the Kronecker product
explicitly. In the product
${\hat{J}_{Q_{k},x}}^{+}\otimes{\hat{J}_{Q_{k},y}}^{-}$ we then use the
permutation equivalence $A\otimes B=PB\otimes AP^{T}$ ($A,B$ square matrices)
to order the projectors $|i\rangle\langle i|_{x}|j\rangle\langle j|_{y}$ to
the left,

$\displaystyle{\hat{J}_{Q_{k},\vec{x}}}^{+}$
$\displaystyle\otimes\hat{J}_{Q_{k},\vec{y}}^{-}=\sum_{i,j}|i\rangle\langle
i|_{x}\otimes J^{(i)+}_{Q_{k},x}\otimes|j\rangle\langle j|_{y}\otimes
J^{(j)-}_{Q_{k},y}$ $\displaystyle=\sum_{i,j}|i\rangle\langle i|_{x}\otimes
P^{T}PJ^{(i)+}_{Q_{k},x}\otimes|j\rangle\langle j|_{y}P^{T}P\otimes
J^{(j)-}_{Q_{k},y}$ $\displaystyle=P^{T}\sum_{i,j}|i\rangle\langle
i|_{x}\otimes|j\rangle\langle j|_{y}\otimes J^{(i)+}_{Q_{k},x}\otimes
J^{(j)-}_{Q_{k},y}P\,.$ (38)

The permutation matrix $P$ then acts on the product state
$\in{\mathbbm{H}_{\mathfrak{h}}}_{x}\otimes{\mathbbm{H}_{\mathfrak{h}}}_{y}$
as
$P(|i\rangle_{x}|\mathfrak{h}_{i}\rangle_{x}|j\rangle_{y}|\mathfrak{h}_{j}\rangle_{y})=|i\rangle_{x}|j\rangle_{y}|\mathfrak{h}_{i}\rangle_{x}|\mathfrak{h}_{j}\rangle_{y}$.
It is in this sense that Hamiltonian Eq.(37) is equivalent to Eq.(18).444We
note that the block-diagonal shape of the local Hamiltonian Eq.(37) may be
interesting to study on its own: In Ref.Gu  _et al._ (2021) such Hamiltonians
are considered for fast-forwarding time evolution, i.e. a time evolution whose
gate complexity is sublinear in $t$.

By reordering the basis, the ancillary register $|i\rangle$ in Eq.(37) can be
chosen to just encode the sectors $|B|=0,1,2$. In Fig.3 we show schematically
the block-diagonal structure of the ladder operators, after reordering and
zero-padding.

Figure 3:  Block-diagonal structure of the ladder operators $J_{Q_{k}}^{+/-}$
after reordering the basis and zero-padding (red dashed line).

As practically we will use the ancillary states as control state, a unary
encoding of the 3 states seems advisable, thus demanding 3 ancillary bits per
site $x$. The actual state at this site will then be encoded using 6 qubits,
hence totalling 9 qubits to represent $|\mathfrak{h}\rangle_{x}$. For the
sector $i=|B|=1$, we can decompose the $J_{Q_{k}}^{(i)+/-}$ even further (c.f.
Fig.3) by noticing that

$\displaystyle J_{Q_{k}}^{(|B|=1)-/+}=|0\rangle\langle 0|\otimes
J_{Q_{k}}^{(B=-1)-/+}+|1\rangle\langle 1|\otimes J_{Q_{k}}^{(B=1)-/+}\,.$ (39)

Because of the $J_{Q_{k}}^{(B=-1)-/+}=J_{Q_{k}}^{(B=1)-/+}$, the qubit count
of 6 “substate” qubits will be reduced by one, ultimately resulting in a lower
gate count.

Figure 4:  The number of partitions $n^{P}$ for different number Pauli strings
$P_{k}$ in our theory, using a sequential approach. Also shown are the results
obtained for a selected set of quantum chemistry Hamiltonians, along with data
taken from Developers (2023) (Table 3) on the same set, for comparison. Figure
5:  The number of entangling CNOT gates and corresponding gate depth for the
circuits implementing the unitaries Eq.(40) in the baryonic sectors $(i,j)$ at
sites $x$ and $y$, respectively.

Finally, to obtain the circuit depth for the time evolution with the
Hamiltonian Eq.(37), we have to diagonalize and optimize each of the factors

$\displaystyle e^{-i\delta tH_{x,y}}$ $\displaystyle=e^{i\frac{\delta
t}{2}\sum_{ij}|i\rangle\langle i|_{x}|j\rangle\langle
j|_{y}\sum_{Q_{k}}\left({\hat{J}_{Q_{k},x}}^{(i)+}{\hat{J}_{Q_{k},y}}^{(j)-}+{\hat{J}_{Q_{k},x}}^{(i),-}{\hat{J}_{Q_{k},y}}^{(j),+}\right)}$
$\displaystyle=\prod_{ij}|i\rangle\langle i|_{x}|j\rangle\langle
j|_{y}e^{i\frac{\delta
t}{2}\sum_{Q_{k}}\left({\hat{J}_{Q_{k},x}}^{(i)+}{\hat{J}_{Q_{k},y}}^{(j)-}+{\hat{J}_{Q_{k},x}}^{(i),-}{\hat{J}_{Q_{k},y}}^{(j),+}\right)}$
$\displaystyle\equiv U^{(0,0)}U^{(1,0)}U^{(0,1)}U^{(1,1)}\,$ (40)

where the individual Trotter factors for different sectors commute
(orthogonality). The $U^{(i,j)}$-diagonalization can be done using the methods
presented in van den Berg and Temme (2020); Murairi  _et al._ (2022). Given
the problem size involving a high dimensional Hamiltonian with 6-qubit-
to-6-qubit coupling as can be seen in Eq.(40), a challenge that arises during
the optimization is the partitioning into commuting sets of the large number
of Pauli strings $P_{k}$ obtained from the expansion

$\displaystyle\sum_{Q_{k}}\left({\hat{J}_{Q_{k},x}}^{(i)+}{\hat{J}_{Q_{k},y}}^{(j)-}+{\hat{J}_{Q_{k},x}}^{(i),-}{\hat{J}_{Q_{k},y}}^{(j),+}\right)=\sum_{l}^{N_{ij}}c_{l}P_{l}\,.$
(41)

With the number of strings $N_{ij}$ being of $\mathcal{O}(10^{6})$ we chose a
sequential approach van den Berg and Temme (2020) over graph coloring. Fig.4
shows the resulting number of partitions $n^{P}_{ij}$ for the different
$N_{ij}$ appearing in our problem, along with a selected set of partitions
obtained from quantum chemistry Hamiltonians available from Developers (2023).
Unsurprisingly, the $n^{P}_{ij}$ do not only depend on the number of Pauli
strings (in the quantum chemistry problems shown the number of qubits varies
between $4-48$), but also on their structure, which often involves identities
$I$ for the molecules shown, thus increasing their “commutativity”. Once the
partitioning is obtained, we proceed with diagonalization and optimization van
den Berg and Temme (2020); Murairi  _et al._ (2022) and ultimately arrive at
the circuits to simulate the sectors of the present theory, along with their
gate count, using IBM Qiskit Qiskit contributors (2023). In Fig.5 we display
the circuit depth and CNOT-count of the circuits implementing the unitary
$U^{(i,j)}$ of Eq.(40), as expected the cost is largely dominated by the
purely mesonic sector with unitary $U^{(0,0)}$ owing to its large dimension.
This statement remains true, even if we include the parameter regime of
$\mu_{B},\mu_{I}>0$: The time evolution involving $\hat{\mathcal{N}}_{B}$ and
$\hat{\mathcal{N}}_{I}$ (defined immediately below Eq.(18)) will again be
controlled by the ancilla register encoding the baryonic sectors $|B|>0$. As
both operators, $\hat{\mathcal{N}}_{B}$ and $\hat{\mathcal{N}}_{I}$, are
diagonal in the chosen basis, they can be readily decomposed into sequences of
CNOT gates and $R_{Z}$-gates of modest length, without Trotterization. The
resource estimate given in Fig.5, along with the qubit requirement of 9 qubits
per site $x$ (3 ancillary qubits for an unary encoding of the baryonic sector
and 6 qubits encoding the actual state) hence represent the central result of
this section.

####  III.2.2 Mapping to qudits

Turning to $d$-level systems with $d>2$, we can once again rewrite Eq.(37) in
a form using the hermitian operators $J^{1/2}$ for each pion, to obtain

$\displaystyle\hat{\mathcal{H}}$ $\displaystyle\sim\sum_{\langle
x,y\rangle}\sum_{i,j}|i\rangle\langle i|_{x}|j\rangle\langle
j|_{y}\sum_{Q_{k}}\left({\hat{J}_{Q_{k,x}}}^{1,(i)}{\hat{J}_{Q_{k,y}}}^{1,(j)}+{\hat{J}_{Q_{k,x}}}^{2,(i)}{\hat{J}_{Q_{k,y}}}^{2,(j)}\right).$
(42)

The labelling of the baryon sectors $|B|=0,1,2$ then suggests a mixed qudit
architecture as was investigated for state-preparation purposes in Mato  _et
al._ (2024), s.t. one can use a qutrit as local ancilla controlling the
baryonic sectors and a qudit with $d=50$ to encode the actual state at site
$x$. To require a qudit of this high dimensionality may seem premature at
first glance. However, recent progress in trapped-ion systems with ^137Ba^+ as
information carrier Low  _et al._ (2023) demonstrated control of a qudit with
up to $d=13$ states, anticipating encodings up to $d=25$ which already
corresponds to half the number of basis states in the $B=0$ sector of our
theory and is capable to completely encode the $B>1$ sectors. Moreover, in
Dong  _et al._ (2023) a photonic qudit with $d=25$ was realized using a cold-
atom approach, showing that such high-dimensional Hilbert spaces can already
be experimentally realized on different platforms. We hence anticipate the
availability of qudit platforms with $d\lesssim 50$ for our purposes and favor
this over the natural possibility of representing the $\mathcal{O}(50)$ states
by the product states $|i_{1}\rangle\otimes|i_{2}\rangle$ of two qudits of
lower dimensionality.

With the above specific architecture in mind, we trotterize the time evolution
using the Hamiltonian Eq.(42) by breaking up the exponential for each pion and
represent the $J^{1,(i)}_{Q_{k}}$ and $J^{2,(i)}_{Q_{k}}$ in terms of the
$\mathcal{X}$ and $\mathcal{Y}$, Eq.(22), in complete analogy to the one-
flavor case,

$\displaystyle
J^{1,(i)}_{Q_{k}}=\sum_{l}c_{l}\mathcal{X}_{m_{l}n_{l}},\,J^{2,(i)}_{Q_{k}}=\sum_{l}d_{l}\mathcal{Y}_{m_{l}n_{l}},\forall\,i,k\,.$
(43)

Depending on the grouping of the $\mathcal{X}$ and $\mathcal{Y}$ into
generalized MS gates, it is conceivable that the entangling gate count is
reduced drastically in each sector combination $(i,j)$ when compared to the
qubit case. After partial resummation of the $\mathcal{X}_{m_{l}n_{l}}$ and
$\mathcal{Y}_{m_{l}n_{l}}$ in Eq.(43) we obtain

$\displaystyle
J^{1,(i)}_{Q_{k}}=\sum^{N_{i}}_{n}\alpha^{(i,k)}_{n}A^{(i,k)}_{n},\,J^{2,(i)}_{Q_{k}}=\sum^{N_{i}}_{n}\beta^{(i,k)}_{n}B^{(i,k)}_{n}\,,$
(44)

where the $A_{n}$ and $B_{n}$ represent linear combinations (c.f. Eq.(23)) of
the $\mathcal{X}$ and $\mathcal{Y}$, respectively. One mesonic Trotter step
would hence involve the unitaries

$\displaystyle e$ ${}^{-i\delta tH_{x,y}}\approx\prod_{i,j}|i\rangle\langle
i|_{x}|j\rangle\langle j|_{y}\times$ $\displaystyle\prod_{k,n,m}$
$\displaystyle e$ ${}^{-i\delta
t\alpha^{(i,k)}_{n}\alpha^{(j,k)}_{m}A^{(i,k)}_{n}\otimes
A^{(j,k)}_{m}}e^{-i\delta
t\beta^{(i,k)}_{n}\beta^{(j,k)}_{m}B^{(i,k)}_{n}\otimes B^{(j,k)}_{m}}\,.$

Considering that, in the $|B|$=1 sector, there are 20 transitions per pion
current $Q_{k}$, this suggests $N_{1}$ = 5 generalized MS gates à 4
transitions. Employing the same principle for $B=0$, yields $N_{0}=7$
generalized MS gates to encode 28 transitions, c.f. Fig.1. In Table 2 we
summarize the arising gate count in terms of generalized MS-gates for fixed
combination of baryonic sectors $(i,j)$ and fixed pion current $Q_{k}$.

Baryon Sector | 0 | 1
---|---|---
0 | $7\times 7\times 2$ | $7\times 5\times 2$
1 | $7\times 5\times 2$ | $5\times 5\times 2$
Table 2:  Estimated count of generalized MS gates for the baryon sector
combination $(i,j)$ per pion current $Q_{k}$. We assumed that the $28$
transitions in the $B=0$ sector can be partitioned into 7 gates à 4
transitions. By analogy for $|B|=1$, we decompose the 20 transitions into 5
gates.

To complete the estimate for time evolving the two-flavor theory with trapped-
ion qudits, we note on the side that the diagonal, local terms
$\hat{\mathcal{N}}_{B}$ and $\hat{\mathcal{N}}_{I}$ in Eq.(17) can be treated
as in the one-flavor case: If we decompose
$\hat{\mathcal{N}}_{B}=\sum_{x,i}|i\rangle\langle i|\otimes\omega_{x,i}$,
where $\omega_{x,0}=0_{50\times 50}$, $\omega_{x,1}=-I_{20\times 20}\oplus
I_{20\times 20}\oplus 0_{10\times 10}=-\sum_{i=0}^{19}\mathcal{Z}_{i,20+i}$
and $\omega_{x,2}=-2\mathcal{Z}_{0,1}$, the gate decomposition in terms of
diagonal single qudit rotations becomes visible, where a similar decomposition
holds for $\hat{\mathcal{N}}_{I}$.

####  III.2.3 Mapping to qumodes

In view of the 92-dimensional state space of the two-flavor theory (see
Sect.II), mapping the theory to qumode degrees of freedom with a theoretically
infinite dimensional Hilbert space appears promising. In fact, one could be
tempted to simply generalize the relation Eq.(31) obtained for $N_{f}=1$ by
defining
$\hat{J}_{Q_{k},x}^{(i)+}=a^{\dagger(i,k)}_{1,x}\hat{a}^{(i,k)}_{2,x}$ to
arrive at

$\displaystyle\mathcal{H}\sim\sum_{\langle
x,y\rangle}\sum_{i,j}|i\rangle\langle i|_{x}|j\rangle\langle
j|_{y}\sum_{Q_{k}}\left(\hat{a}^{\dagger(i,k)}_{1,x}\hat{a}^{(i,k)}_{2,x}\hat{a}^{\dagger(j,k)}_{2,y}\hat{a}^{(j,k)}_{1,y}+h.c.\right).$

However, this expression does not reproduce the correct matrix elements.
Contrary to the one-flavor case, where the $J^{+/-}$ simply fulfilled the
algebra of a spin-$N_{c}/2$ representation (c.f. Sect.II), the
$\hat{J}_{Q_{k},x}^{(i)+}$ in the $N_{f}=2$-theory are reducible and decompose
into product representations whose structure still has to be determined Unger
and Pattanaik (2022a, b). It is conceivable that the Jordan-Schwinger-map for
the theory with $N_{f}=2$ requires several qumode operators $a_{i}^{k}$ per
pion current $Q_{k}$ and – following the discussion in Sect.III.1.3 – is
likely to be expensive in terms of the number of gates which concludes our
discussion for the two-flavor case.

##  IV Conclusion and Outlook

Table 3:  Storage requirement for lattice volume $N$

Information Carrier |  $N_{f}=1$ | $N_{f}=2$
---|---|---
Qubit |  $3N$ |  $(3+6)N$*
Qudit $(d>2)$ |  $N$ |  $N+N$**
Qumode |  $2N+N$** | $>2N$

  * *

We used unary encoding of the three ancilla states representing the baryon
sectors $|B|$.

  * **

With mixed architecture, i.e. either mixed-qudit or qubit-qumode.

Table 4: Entangling gate count for one mesonic $nn$-Trotter step

Information Carrier |  $N_{f}=1$ | $N_{f}=2$
---|---|---
Qubit |  $\mathcal{O}(10)$ | $\mathcal{O}(10^{6})$
Qudit $(d>2)$ | 2 | $\mathcal{O}(10^{2})$
Qumode |  $\mathcal{O}(10^{2})-\mathcal{O}(10^{3})$* | –

  * *

Depending on the availability of the quartic gate.

We summarize the relevant information on the resource requirement for the
digitization of strong coupling lattice QCD with $N_{f}=1$ and $N_{f}=2$
flavors of staggered quarks in Table 3 and Table 4. As qudits with $d>2$ offer
the possibility to represent a multidimensional (local) state space, their
resource requirements in terms of storage are favorable when compared to
2-level systems and qumodes. For the latter it should be noted that the
scaling of course depends on the particular mapping chosen to encode the
local, discrete state space of our models (angular momentum states) in the
Fock-states of quantum harmonic oscillators, where we used the Jordan-
Schwinger-map. For continuous variables, a different encoding may result in a
more efficient scaling, as is the case for the non-linear sigma model Jha  _et
al._ (2023a, b). The entangling gate count in the qumode case (Table 4) shows
that the availability of non-Gaussian gates (in particular the quartic gate)
has a high impact on the scaling. It remains to be seen if a measurement based
approach as advocated in Abel  _et al._ (2024) offers a practical alternative.
Conversely, for the qudit-case, it was the simultaneous driving of transitions
via generalized MS-gates that allowed to reduce the gate count by at least one
order of magnitude – a reduction that is achieved by more complex experimental
control. With qubits as information carrier, one main driver of gate count is
the partitioning step for large Hamiltonians. Generally graph coloring is
preferred which quickly becomes too expensive with the library we used Hagberg
_et al._ (2008) and we thus reverted to sequential approaches van den Berg and
Temme (2020). This partitioning or clustering is an active field of research
and progress driven by the different optimization goals (e.g. simultaneous
measurement or gate count reduction as in our case) will likely reduce the
qubit result for $N_{f}=2$ in Table 4.

As our study solely considers the resource requirements from an implementation
point of view, an extension could take into account the influence of errors on
different platforms and existing error mitigation techniques. Another aspect
that certainly deserves attention is the question of state preparation on
different devices. In particular, for the photonic case, preparing states with
different photon occupation number to encode the mesonic or baryonic sector,
is an interesting subject of study.

###### Acknowledgements.

The work is supported by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) through the grant CRC-TR 211 “Strong-interaction matter
under extreme conditions” – project number 315477589 – TRR 211 and by the
State of Hesse within the Research Cluster ELEMENTS (Project ID 500/10.006).
M.F. acknowledges support by the Munich Institute for Astro-, Particle and
BioPhysics (MIAPbP) which is funded by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) under Germany’s Excellence Strategy –
EXC-2094 – 390783311.

## References

  * Lloyd (1996) S. Lloyd, Science 273, 1073 (1996).
  * Bauer  _et al._ (2022) C. W. Bauer, Z. Davoudi, A. B. Balantekin, T. Bhattacharya, M. Carena, W. A. de Jong, P. Draper, A. El-Khadra, N. Gemelke, M. Hanada, D. Kharzeev, H. Lamm, Y.-Y. Li, J. Liu, M. Lukin, Y. Meurice, C. Monroe, B. Nachman, G. Pagano, J. Preskill, E. Rinaldi, A. Roggero, D. I. Santiago, M. J. Savage, I. Siddiqi, G. Siopsis, D. V. Zanten, N. Wiebe, Y. Yamauchi, K. Yeter-Aydeniz, and S. Zorzetti, (2022), arXiv:2204.03381 [quant-ph] .
  * Di Meglio  _et al._ (2023) A. Di Meglio, K. Jansen, I. Tavernelli, C. Alexandrou, S. Arunachalam, C. W. Bauer, K. Borras, S. Carrazza, A. Crippa, V. Croft, R. de Putter, A. Delgado, V. Dunjko, D. J. Egger, E. Fernandez-Combarro, E. Fuchs, L. Funcke, D. Gonzalez-Cuadra, M. Grossi, J. C. Halimeh, Z. Holmes, S. Kuhn, D. Lacroix, R. Lewis, D. Lucchesi, M. L. Martinez, F. Meloni, A. Mezzacapo, S. Montangero, L. Nagano, V. Radescu, E. R. Ortega, A. Roggero, J. Schuhmacher, J. Seixas, P. Silvi, P. Spentzouris, F. Tacchino, K. Temme, K. Terashi, J. Tura, C. Tuysuz, S. Vallecorsa, U.-J. Wiese, S. Yoo, and J. Zhang, “Quantum computing for high-energy physics: State of the art and challenges. summary of the qc4hep working group,”  (2023).
  * Wilson (1974) K. G. Wilson, Physical Review D 10, 2445 (1974).
  * Kogut and Susskind (1975) J. Kogut and L. Susskind, Phys. Rev. D 11, 395 (1975).
  * Byrnes and Yamamoto (2005) T. Byrnes and Y. Yamamoto, Phys.Rev. A73 (2006) 022328 (2005), 10.1103/PhysRevA.73.022328, arXiv:quant-ph/0510027 [quant-ph] .
  * Tong  _et al._ (2021) Y. Tong, V. V. Albert, J. R. McClean, J. Preskill, and Y. Su, Quantum 6, 816 (2022) 6, 816 (2021), arXiv:2110.06942 [quant-ph] .
  * Davoudi  _et al._ (2020) Z. Davoudi, I. Raychowdhury, and A. Shaw, Phys. Rev. D 104, 074505 (2021) (2020), 10.1103/PhysRevD.104.074505, arXiv:2009.11802 [hep-lat] .
  * Assi and Lamm (2024) B. Assi and H. Lamm,  (2024), 10.48550/ARXIV.2405.12204, arXiv:2405.12204 [hep-lat] .
  * Davoudi  _et al._ (2023) Z. Davoudi, A. F. Shaw, and J. R. Stryker, Quantum 7, 1213 (2023).
  * Wang  _et al._ (2020) Y. Wang, Z. Hu, B. C. Sanders, and S. Kais, Frontiers in Physics 8, 479 (2020), arXiv:2008.00959 [quant-ph] .
  * Weedbrook  _et al._ (2012) C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Reviews of Modern Physics 84, 621 (2012).
  * Lloyd and Braunstein (1999) S. Lloyd and S. L. Braunstein, Physical Review Letters 82, 1784 (1999).
  * Fromm  _et al._ (2023) M. Fromm, O. Philipsen, W. Unger, and C. Winterowd, EPJ Quantum Technology 11 (2023), 10.1140/epjqt/s40507-024-00236-y, arXiv:2308.03196 [hep-lat] .
  * Rossi and Wolff (1984) P. Rossi and U. Wolff, Nuclear Physics B 248, 105 (1984).
  * Karsch and Mütter (1989) F. Karsch and K.-H. Mütter, Nuclear Physics B 313, 541 (1989).
  * de Forcrand and Fromm (2010) P. de Forcrand and M. Fromm, Physical Review Letters 104, 112005 (2010).
  * de Forcrand  _et al._ (2014) P. de Forcrand, J. Langelage, O. Philipsen, and W. Unger, Physical Review Letters 113, 152002 (2014).
  * Klegrewe and Unger (2020) M. Klegrewe and W. Unger, Physical Review D 102, 034505 (2020).
  * Unger and Pattanaik (2022a) W. Unger and P. Pattanaik, in  _Proceedings of The 38th International Symposium on Lattice Field Theory — PoS(LATTICE2021)_, LATTICE2021 (Sissa Medialab, 2022).
  * Unger and Pattanaik (2022b) W. Unger and P. Pattanaik, in  _Proceedings of The 39th International Symposium on Lattice Field Theory — PoS(LATTICE2022)_, LATTICE2022 (Sissa Medialab, 2022).
  * Smit (1980) J. Smit, Nuclear Physics B 175, 307 (1980).
  * Unger (2024) W. Unger, in  _Proceedings of The 40th International Symposium on Lattice Field Theory — PoS(LATTICE2023)_, LATTICE2023 (Sissa Medialab, 2024).
  * O’Leary  _et al._ (2006) D. P. O’Leary, G. K. Brennen, and S. S. Bullock, Physical Review A 74, 032334 (2006).
  * Molmer and Sorensen (1998) K. Molmer and A. Sorensen, Phys.Rev.Lett. 82 (1999) 1835-1838 82, 1835 (1998), arXiv:quant-ph/9810040 [quant-ph] .
  * Ringbauer  _et al._ (2021) M. Ringbauer, M. Meth, L. Postler, R. Stricker, R. Blatt, P. Schindler, and T. Monz, Nature Physics 18, 1053 (2022) 18, 1053 (2021), arXiv:2109.06903 [quant-ph] .
  * Low  _et al._ (2019) P. J. Low, B. M. White, A. A. Cox, M. L. Day, and C. Senko, Phys. Rev. Research 2, 033128 (2020) 2, 033128 (2019), arXiv:1907.08569 [physics.atom-ph] .
  * Calajò  _et al._ (2024) G. Calajò, G. Magnifico, C. Edmunds, M. Ringbauer, S. Montangero, and P. Silvi,  (2024), 10.48550/ARXIV.2402.07987, arXiv:2402.07987 [quant-ph] .
  * Illa  _et al._ (2024) M. Illa, C. E. P. Robin, and M. J. Savage,  (2024), 10.48550/ARXIV.2403.14537, arXiv:2403.14537 [quant-ph] .
  * Schwinger (1952) J. Schwinger,  _On Angular Momentum_ , Tech. Rep. (1952) report NYO-3071.
  * Annabestani  _et al._ (2020) R. Annabestani, B. Gupt, and B. R. Bardhan,  (2020), 10.48550/ARXIV.2009.09455, arXiv:2009.09455 [quant-ph] .
  * Blais  _et al._ (2020) A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, Rev. Mod. Phys. 93, 25005 (2021) 93, 025005 (2020), arXiv:2005.12667 [quant-ph] .
  * Killoran  _et al._ (2019) N. Killoran, J. Izaac, N. Quesada, V. Bergholm, M. Amy, and C. Weedbrook, Quantum 3, 129 (2019).
  * Gu  _et al._ (2021) S. Gu, R. D. Somma, and B. Şahinoğlu, Quantum 5, 577 (2021).
  * Developers (2023) T. Q. N. Developers, “Qiskit nature 0.6.0,”  (2023).
  * van den Berg and Temme (2020) E. van den Berg and K. Temme, Quantum 4, 322 (2020) (2020), 10.22331/q-2020-09-12-322, arXiv:2003.13599 [quant-ph] .
  * Murairi  _et al._ (2022) E. M. Murairi, M. J. Cervia, H. Kumar, P. F. Bedaque, and A. Alexandru, (2022), arXiv:2208.11789 [hep-lat] .
  * Qiskit contributors (2023) Qiskit contributors, “Qiskit: An open-source framework for quantum computing,”  (2023).
  * Mato  _et al._ (2024) K. Mato, S. Hillmich, and R. Wille, Design Automation Conference 2024, San Francisco (USA) (2024), 10.1145/3649329.3656260, arXiv:2406.03531 [quant-ph] .
  * Low  _et al._ (2023) P. J. Low, B. White, and C. Senko,  (2023), 10.48550/ARXIV.2306.03340, arXiv:2306.03340 [quant-ph] .
  * Dong  _et al._ (2023) M.-X. Dong, W.-H. Zhang, L. Zeng, Y.-H. Ye, D.-C. Li, G.-C. Guo, D.-S. Ding, and B.-S. Shi, Physical Review Letters 131, 240801 (2023).
  * Jha  _et al._ (2023a) R. G. Jha, F. Ringer, G. Siopsis, and S. Thompson,  _Proceedings of The 40th International Symposium on Lattice Field Theory — PoS(LATTICE2023)_ ,  LATTICE2023 (2023a), 10.22323/1.453.0230, arXiv:2308.06946 [hep-lat] .
  * Jha  _et al._ (2023b) R. G. Jha, F. Ringer, G. Siopsis, and S. Thompson, Physical Review A 109, 052412 (2023b), arXiv:2310.12512 [quant-ph] .
  * Abel  _et al._ (2024) S. Abel, M. Spannowsky, and S. Williams,  (2024), 10.48550/ARXIV.2403.10619, arXiv:2403.10619 [quant-ph] .
  * Hagberg  _et al._ (2008) A. A. Hagberg, D. A. Schult, and P. J. Swart, in  _Proceedings of the 7th Python in Science Conference_ , edited by G. Varoquaux, T. Vaught, and J. Millman (Pasadena, CA USA, 2008) pp. 11 – 15\.
  * Kalajdzievski and Arrazola (2018) T. Kalajdzievski and J. M. Arrazola, Phys. Rev. A 99, 022341 (2019) 99, 022341 (2018), arXiv:1811.10651 [quant-ph] .

##  Appendix A Digitization using qumodes

###  A.1 Exact gate decomposition using qumodes, gate count

All of the eight terms appearing in Eq.(32) can be rewritten as
$\hat{x}_{1,x}\hat{x}_{2,x}\hat{x}_{1,y}\hat{x}_{2,y}$ by means of unitary
conjugation with Fourier gates, s.t. for example

$\displaystyle e^{i\delta
t\hat{x}_{1,x}\hat{p}_{2,x}\hat{x}_{1,y}\hat{x}_{2,y}}=F_{2,x}e^{i\delta
t\hat{x}_{1,x}\hat{x}_{2,x}\hat{x}_{1,y}\hat{x}_{2,y}}F^{\dagger}_{2,x}\,.$
(45)

Next we employ the identity

$\displaystyle\hat{x}_{1,x}\hat{x}_{2,x}\hat{x}_{1,y}\hat{x}_{2,y}=\frac{1}{192}$
$\displaystyle[$
$\displaystyle(x_{1,x}+\hat{x}_{2,x}+\hat{x}_{1,y}+\hat{x}_{2,y})^{4}+(x_{1,x}+\hat{x}_{2,x}-\hat{x}_{1,y}-\hat{x}_{2,y})^{4}$
(46) $\displaystyle+$
$\displaystyle(x_{1,x}-\hat{x}_{2,x}+\hat{x}_{1,y}-\hat{x}_{2,y})^{4}+(x_{1,x}-\hat{x}_{2,x}-\hat{x}_{1,y}+\hat{x}_{2,y})^{4}$
$\displaystyle-$
$\displaystyle(x_{1,x}-\hat{x}_{2,x}-\hat{x}_{1,y}-\hat{x}_{2,y})^{4}-(x_{1,x}+\hat{x}_{2,x}+\hat{x}_{1,y}-\hat{x}_{2,y})^{4}$
$\displaystyle-$
$\displaystyle(x_{1,x}+\hat{x}_{2,x}-\hat{x}_{1,y}+\hat{x}_{2,y})^{4}-(x_{1,x}-\hat{x}_{2,x}+\hat{x}_{1,y}+\hat{x}_{2,y})^{4}]\,.$

Each term appearing on the right can be decomposed via unitary conjugation,
s.t.

$\displaystyle e^{i\delta
t(\hat{x}_{1,x}+\hat{x}_{2,x}+\hat{x}_{1,y}+\hat{x}_{2,y})^{4}}=e^{2i\hat{p}_{1,x}\hat{x}_{2,y}}e^{2i\hat{p}_{1,x}\hat{x}_{1,y}}e^{2i\hat{p}_{1,x}\hat{x}_{2,x}}e^{i\delta
t\hat{x}^{4}_{1,x}}e^{-2i\hat{p}_{1,x}\hat{x}_{2,x}}e^{-2i\hat{p}_{1,x}\hat{x}_{1,y}}e^{-2i\hat{p}_{1,x}\hat{x}_{2,y}}\,,$
(47)

which can be seen by making repeated use of the Baker-Campbell-Hausdorff
formula $e^{A}Be^{-A}=B+[A,B]+\frac{1}{2}[A,[A,B]]+\ldots$. We note that
$e^{i\hat{p}_{k}\hat{x}_{l}}=F_{k}e^{i\hat{x}_{k}\hat{x}_{l}}F^{\dagger}_{k}$,
where the operations on the $rhs$ are already part of the instruction set. The
hopping term Eq.(45) is hence already decomposed into elementary operations,
if the quartic gate $e^{i\delta t\hat{x}^{4}_{j}}$, appearing in Eq.(47), is
part of the instruction set. Otherwise the latter needs to decomposed by means
of an ancilla qumode,
$\hat{x}_{j}^{4}=(\hat{x}_{j}^{2}+\hat{x}_{a})^{2}-\hat{x}_{a}^{2}-2\hat{x}_{a}\hat{x}_{j}^{2}$,
where we can further use that

$\displaystyle e^{i\delta
t(\hat{x}_{j}^{2}+\hat{x}_{a})^{2}}=e^{2i\hat{p}_{a}\hat{x}_{j}^{2}}e^{i\delta
t\hat{x}_{a}^{2}}e^{-2i\hat{p}_{a}\hat{x}_{j}^{2}}$ (48)

and employ the identity Kalajdzievski and Arrazola (2018)

$\displaystyle
e^{3i\alpha^{2}t\hat{p}_{a}\hat{x}_{j}^{2}}=e^{2i\alpha\hat{x}_{j}\hat{x}_{a}}e^{it\hat{p}_{a}^{3}}e^{-i\alpha\hat{x}_{j}\hat{x}_{a}}e^{-it\hat{p}_{a}^{3}}e^{-2i\alpha\hat{x}_{j}\hat{x}_{a}}e^{it\hat{p}_{a}^{3}}e^{i\alpha\hat{x}_{j}\hat{x}_{a}}e^{-it\hat{p}_{a}^{3}}e^{i\frac{3}{4}\alpha^{3}t\hat{x}_{j}^{3}},\alpha\in\mathbb{R}\,.$
(49)

A decomposition of the quartic gate $Q(\delta t)=e^{i\delta t\hat{x}^{4}_{j}}$
using the relations Eqs.(48) and (49) turns out to be costly with a total of
55 elementary operations, of both, single qumode and entangling two-qumode
type. This will ultimately drive the main cost of a mesonic time evolution.
Consequently, in addition to the quartic gate $Q$, the decomposition shown in
Eq.(47) will take 6 $CX$-gates which will amount to 4 Fourier gates $F$ and 6
$CZ$-gates of our instruction set Eq.(33). As each term appearing on the $rhs$
of Eq.(46) can be decomposed analogously, its exponentiation will lead to
massive cancellations and combinations, resulting in a total count of
$8Q,19CZ$ and $18F$ for the exponential of Eq.(46).
We now trotterize the exponential of the mesonic hopping term in Eq.(32) into
the eight non-commuting terms of the form Eq.(45) with appropriate
combinations of Fourier gates $p(F)$, i.e.

$\displaystyle e^{i\delta
t\sum_{\\{i_{k}\\}}(\hat{q}^{i_{1}}_{1,x}\hat{q}^{i_{2}}_{2,x}\hat{q}^{i_{3}}_{2,y}\hat{q}^{i_{4}}_{1,y})s_{\\{i_{1}\ldots
i_{4}\\}}}=\prod_{\\{i_{k}\\}}p(F)e^{i\delta
t\hat{x}_{1,x}\hat{x}_{2,x}\hat{x}_{1,y}\hat{x}_{2,y}s_{\\{i_{1}\ldots
i_{4}\\}}}p(F^{\dagger})+\mathcal{O}(\delta t^{2})\,.$ (50)

Given our considerations above, this can be done with cancellations in the
appearing Fourier gates and CZ gates, s.t. the total gate count for the above
operation is given by $64Q,147CZ,149F$. If we consider the single qubit gate Q
as one operation, then the corresponding gate depth is given by 348, as only
interjacent Fourier gates can be carried out in parallel.