Discrete Cylindrical Vector Beam Generation from an Array of Optical Fibers

R. Steven Kurti1,∗, Klaus Halterman1, Ramesh K. Shori,1 and Michael J.
Wardlaw2

1Research Division, Physics and Computational Sciences Branch, Naval Air
Warfare Center,
China Lake, California, 93555, USA
2Office of Naval Research, Arlington, Virgina, 22203, USA.
∗Corresponding author<EMAIL_ADDRESS>
###### Abstract

A novel method is presented for the beam shaping of far field intensity
distributions of coherently combined fiber arrays. The fibers are arranged
uniformly on the perimeter of a circle, and the linearly polarized beams of
equal shape are superimposed such that the far field pattern represents an
effective radially polarized vector beam, or discrete cylindrical vector (DCV)
beam. The DCV beam is produced by three or more beams that each individually
have a varying polarization vector. The beams are appropriately distributed in
the near field such that the far field intensity distribution has a central
null. This result is in contrast to the situation of parallel linearly
polarized beams, where the intensity peaks on axis.

OCIS codes: 140.3298, 140.3300, 060.3510.

The propagation of electromagnetic fields from multiple fibers can result in
complex far-field intensity profiles that depend crucially on the individual
near field polarization and phase. Coherently phased arrays have been used in
defense and communications systems for many years, where the antenna ensemble
forms a diffraction pattern that can be altered by changing the antenna
spacing, amplitude, and relative phase relationships. In the optical regime,
similar systems have recently been formed by actively or passively locking the
phases of two or more identical optical beams [1, 2, 3, 4]. While coherently
combined fiber lasers are increasingly gaining acceptance as sources for high
power applications, most approaches involve phasing a rectangular or hexagonal
grid of fibers with collimated beams that are linearly polarized along the
same axis [5]. The resultant far field diffraction pattern therefore is peaked
in the center. If the polarization is allowed to vary, as it does in radial
vector beams, the diffraction pattern vanishes in the center[6], which can
result in a myriad of device applications, discussed below.

Several types of polarization states have been investigated, including radial
and azimuthal [7, 8, 6, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18]. These beams
are used in mitigating thermal effects in high power lasers [19, 20, 21],
laser machining [22, 23], and particle acceleration interactions [24, 25].
They can even be used to generate longitudinal electric fields when tightly
focused [26, 27], and although they are typically formed in free space laser
cavities using a conical lens, or axicon, they can also be formed and guided
in fibers [28, 29, 30]. Typical methods for creating radially polarized beams
fall into two categories: the first involves an intracavity axicon in a laser
resonator to generate the laser mode, the second begins with a single beam and
rotates the polarization of portions of the beam to create an inhomogeneously
(typically radially or azimuthally) polarized beam. In this Letter we
demonstrate a novel approach that takes an array of Gaussian beams, each with
the appropriately oriented linear polarization, and then superimpose them in
the far-field, effectively creating a DCV beam.

Our starting point in determining the propagation of polarized electromagnetic
beams is the vector Helmholtz equation, $\nabla\times\nabla\times{\bf
E}({\boldsymbol{\rho}},z)-k_{0}^{2}{\bf E}({\boldsymbol{\rho}},z)=0$. We use
cylindrical coordinates, so that $\boldsymbol{\rho}$ is the usual transverse
coordinate, $k_{0}=\omega/c$, and the usual sinusoidal time dependence has
been factored out. Applying Lagrange’s formula permits the wave equation to be
written in terms of the vector Laplacian: $\nabla^{2}{\bf
E}({\boldsymbol{\rho}},z)+k_{0}^{2}{\bf E}({\boldsymbol{\rho}},z)=0$, which
can then be reduced to the paraxial wave equation. The paraxial solutions are
then inserted into the familiar Fraunhofer diffraction integral, which for
propagation at sufficiently large $z$, yields the electric field,

$\displaystyle{\bf
E}(k_{x},k_{y})=\exp\Bigl{(}\frac{ik_{0}\rho^{2}}{2z}\Bigr{)}\frac{\exp(ik_{0}z)}{i\lambda
z}$
$\displaystyle\times\sum_{j=1}^{N}\int_{0}^{a}rdr\int_{0}^{2\pi}d\phi\,e^{-i{\bf
k}\cdot({\bf r}+{\bf s}_{j})}{\bf E}^{0}({\bf r}+{\bf s}_{j},0),$ (1)

where ${\bf E}^{0}({\bf r}+{\bf s}_{j},0)$ is the incident electric field in
the plane of the fiber array.

Fig. 1: Example of the multiple fiber setup. A specific case of $N=3$ is
shown, where each beam is represented by a circle of radius $a$ and arranged
uniformly on a circle of radius $R$. Points within the $i$th hole are located
by ${\bf d}={\bf s}_{i}+{\bf r}$, with the vector ${\bf r}$ originating at the
center of each hole.

It is clear that longitudinal polarization is not considered here, consistent
with the paraxial approximation. The far field intensity will be a complex
diffraction pattern that depends on the individual beam intensity profile in
the near field as well its polarization state and optical phase. A diagram
illustrating the geometry used is shown in Fig. 1, where an array of three
circular holes of radius $a$ are equally distributed on the circumference of a
circle of radius $R$. The center of each hole is located at ${\bf
s}_{j}=R\hat{\bf r}_{j}$, where $\hat{\bf r}_{j}\equiv(\hat{\bf
x}\cos\theta_{j}+\hat{\bf y}\sin\theta_{j})$ is the unit radial vector, and
${\bf r}$ is the coordinate relative to each center. The center of each fiber
is separated by an angle $\theta_{j}=2\pi(j-1)/N$. One can scale to any number
of beams, constrained only by the radius $R$. For identical Gaussian beams,
the incident field is expressed as, ${\bf E}^{0}={E}^{0}({\bf r})\hat{\bf
r}_{j}$, where ${E}^{0}({\bf r})=E_{0}e^{-r^{2}/w_{0}^{2}}$. This permits Eq.
(S0.Ex1) to be separated,

$\displaystyle{\bf E}(k_{x},k_{y})={\cal F}(\rho,z)\sum_{j=1}^{N}e^{-i{\bf
k}\cdot{\bf s}_{j}}\hat{\bf r}_{j},$ (2)

where ${\cal
F}(\rho,z)=E_{0}k_{0}/z\int^{a}_{0}rdr\,J_{0}(k_{0}r\rho/z)e^{-r^{2}/w_{0}^{2}}$,
and the prefactors that do not contribute to the intensity have been
suppressed. It is clear from Eq. (2) that the far field transform preserves
the radially symmetric polarization state of the system, as expected. Note
that in the limit $w_{0}/a\gg 1$, and using the relation
$xJ_{0}(x)=[xJ_{1}(x)]^{\prime}$, we have ${\cal F}(\rho,z)\approx
E_{0}(a/\rho)J_{1}(k_{0}a\rho/z)$, which gives the expected Fraunhofer
diffraction pattern for a circular aperture of radius $a$. In the opposite
limit, where the Gaussian profile is narrow enough so that the aperture
geometry has little effect, we have ${\cal F}(\rho,z)\approx
E_{0}k_{0}w_{0}^{2}/(2z)\exp[-(k_{0}\rho w_{0}/(2z))^{2}]$. For the multiple
beam arrangements investigated, corresponding to $N=3,4,6$, we then can
calculate the intensity, $I_{N}$, explicitly:

$\displaystyle{\cal K}_{3}(\phi)$
$\displaystyle=3-\cos({\sqrt{3}}k_{y}R)-2\cos(\frac{\sqrt{3}}{2}k_{y}R)\cos\Bigl{(}\frac{3}{2}k_{x}R\Bigr{)},$
(3) $\displaystyle{\cal K}_{4}(\phi)$
$\displaystyle=4\bigl{[}\sin^{2}(k_{x}R)+\sin^{2}(k_{y}R)\bigr{]},$ (4)
$\displaystyle{\cal K}_{6}(\phi)$
$\displaystyle=4\Bigl{[}\Bigl{(}\cos\Bigl{(}\frac{\sqrt{3}}{2}k_{y}R\Bigr{)}\sin\Bigl{(}\frac{1}{2}k_{x}R\Bigr{)}+\sin\bigl{(}k_{x}R\bigr{)}\Bigr{)}^{2}$
$\displaystyle\,\,\,\,\,\,+3\cos^{2}\Bigl{(}\frac{1}{2}k_{x}R\Bigr{)}\sin^{2}\Bigl{(}\frac{\sqrt{3}}{2}k_{y}R\Bigr{)}\Bigr{]},$
(5)

where $I_{N}\equiv{\cal F}^{2}(\rho,z){\cal K}_{N}(\phi;\rho,z)$,
$k_{x}=(k_{0}/z)\rho\cos\phi$ and $k_{y}=(k_{0}/z)\rho\sin\phi$. To rotate the
array, one can perform a standard rotation ${\bf r}^{\prime}={\cal
R}(\phi^{\prime}){\bf r}$, so e.g., a $\pi/4$ rotation would give, ${\cal
K}_{4}(\phi^{\prime})=4(1-\cos(\sqrt{2}k_{x}R)\cos(\sqrt{2}k_{y}R))$.

We employ two different methods to create the DCV beams previously described.
In each method, the experimental configuration utilizes collimated Gaussian
beams that are in phase. The arrays are cylindrically symmetric both with
respect to intensity and polarization. We focus the beams with a transform
lens in order to simulate the far field intensity at the lens focus. The focal
spot is then imaged onto a camera by a microscope objective in order to fill
the CCD array. The laser source is a $100$ mW continuous wave (CW) Nd:YAG
operating at $\lambda=1.064\mu$m.

Fig. 2: Diagram of the experimental setup for DCV beam generation. The
oscillator output is split and phase modulated (S/PM). The output of each
fiber is collimated by means of a lens, L1, and then the polarization is
rotated by a half-wave plate (HWP). The beam ensemble is then focused to a
point which is imaged by a microscope objective (MO) onto a charge-coupled-
device (CCD) array.

The first method (for the case of $N=4$) involves expanding a beam from the
Nd:YAG by means of a telescope beam expander to a diameter of about 1cm. The
diffractive optical element then creates an $8\times 8$ array of beams. Four
of the beams are reflected and made nearly parallel by a segmented mirror, and
then linearly polarized upon passing through a set of four half-wave plates.
The phases of each beam are made identical by having them pass through
articulated glass slides.

|
---|---
|
|
Fig. 3: Far field intensity profiles of effectively radially polarized beams.
The left panels correspond to the experimental images and the right panels
represent the theoretical results. We consider a 3 (top row), 4 (middle row),
and 6 (bottom row) beam arrangement. For $N=3$, we take $w_{0}=32$ $\mu$m and
$R=5w_{0}$. For $N=4$, $w_{0}=42$ $\mu$m, and $R=2.7w_{0}$. For $N=6$, we have
$w_{0}=43$ $\mu$m and $R=3.8w_{0}$. Note that in general, a radially polarized
beam undergoes a discontinuity at the origin, and thus the intensity must
vanish there.

The second method (for the case of $N=3$ and $6$), which is simpler for larger
arrays, is shown in Fig. 2. In this method the Nd:YAG is coupled into a
polarization maintaining fiber. This signal is then split into multiple copies
by means of a lithium niobate waveguide 8-way splitter and 8-channel electro-
optic modulator. Each path can then have a separate phase modulation to ensure
proper phasing of the beams. These signals are then propagated through fiber
and coupled out and collimated by a 1 in. lens. The lens is slightly
overfilled so the Gaussian outputs of the fibers are truncated at $1.1\times$
(the $e^{-1}$ radius).

Experimental results for the radial vector beam generation found good
agreement with theoretical calculations, as illustrated in Fig. 3. Clearly,
the interference patterns and symmetry correlate well with the corresponding
calculations in the three beam case. In particular, for $N=4$ (middle row),
the central null surrounded by a checkerboard peak structure. This pattern
differs from the linear polarization by a rotation of $\pi/4$, based on a
simple phase argument. The $N=6$ case is shown in the bottom row, where
clearly the intensity vanishes at the origin, followed by the formation of
bright hexagonal rings. The three beam configuration is shown in the top row,
where again there is satisfactory agreement between the calculated and
measured results. Any observed discrepancies may be reduced with an
appropriately incorporated feedback system. A substantial fraction of the
energy in the plots of Fig. 3 is contained within the first regions of high
intensity peaks. We illustrate this for the $N=4$ case where the net intensity
within a circular region of radius $\rho$, $U(\rho)$, is calculated to be
approximately $U(\rho)\approx k_{0}\rho z/(2R)[k_{0}\rho R/z-J_{1}(2k_{0}\rho
R/z)]$. Inserting the appropriate parameters and normalizing to the intensity
integrated over the entire image plane, yields $50\%$ of the distribution is
contained in the neighborhood of the first main peaks.

In conclusion, it has been shown by employing two different configurations,
and by carefully controlling the polarization, a central null can be created
in the far field. In previous work, the central portion of a phased array of
lasers typically has been a peaked function. Now for the first time to our
knowledge, a central null has been formed. Similar results have been shown
from single aperture lasers using an axicon to generate an annular mode inside
a laser cavity. In the present case however, there is direct control over
which mode will be generated from the same system. As for scaling, there is no
fundamental limitation on the number of beams used. Several methods are
available to incorporate a greater number of beams, including concentric beam
placement (following the same beam-combining prescription described in this
Letter) that could scale in roughly the same pattern as a Bessel beam. It is
also possible to use fiber lasers with active phasing, thus potentially
scaling to many hundreds of beams[5]. Our system also has a practical
advantage over typical high power Gaussian beam applications which have the
drawback of having the beam concentrated near the center, where the reflecting
beam director (such as a telescope obscuration) resides. Our device on the
other hand, has no power propagating along the center axis of the beam, and
thus for high power applications requiring a center obscuration telescope beam
director, our proposed system offers new advances.

We thank S. Feng for many useful discussions.

## References

  * [1] “Theory of electronically phased coherent beam combination without a reference beam”, T. M. Shay, Opt. Express 14, 12188 (2006).
  * [2] “Self-organized coherence in fiber laser arrays”, H. Bruesselbach, D. C. Jones, M. S. Mangir, M. I. Minden, and J. L. Rogers, Opt. Lett. 30, 1339, (2003).
  * [3] “Coherent combining of spectrally broadened fiber lasers,” T. B. Simpson, F. Doft, P. R. Peterson, and A. Gavrielides, Opt. Express 15, 11731 (2007)
  * [4] “Coherent addition of fiber lasers by use of a fiber coupler,” Shirakawa, T. Saitou, T. Sekiguchi, and K. Ueda, Opt. Express 10, 1167 (2002).
  * [5] “Self-Synchronous and Self-Referenced Coherent Beam Combination for Large Optical Arrays,” T.M. Shay, V. Benham, J.T. Baker, A.D. Sanchez, D. Pilkington, C.A. Lu, IEEE J. Sel. Top. Quant. El. 13, 480 (2007).
  * [6] “The formation of laser beams with pure azimuthal or radial polarization,” R. Oron, S. Blit, N. Davidson, A. A. Friesem, Z. Bomzon, and E. Hasman, Appl. Phys. Lett. 77, 3322 (2000).
  * [7] “Efficient extracavity generation of radially and azimuthally polarized beams,” G. Machavariani, Y. Lumer, I. Moshe, A. Meir, and S. Jackel, Opt. Lett. 32, 1468 (2007).
  * [8] “Intracavity generation of radially polarized CO2 laser beams based on a simple binary dielectric diffraction grating,” T. Moser, J. Balmer, D. Delbeke, P. Muys, S. Verstuyft, and R. Baets, Appl. Opt. 45, 8517 (2006).
  * [9] “Simple interferometric technique for generation of a radially polarized light beam,” N. Passilly, R. de Saint Denis, K. A t-Ameur, F. Treussart, R. Hierle, and J. -F. Roch, J. Opt. Soc. Am. A 22, 984 (2005).
  * [10] “Generating radially polarized beams interferometrically,” S. C. Tidwell, D. H. Ford, and W. D. Kimura, Appl. Opt. 29, 2234 (1990).
  * [11] “Laser beams with axially symmetric polarization,” A V Nesterov et al., J. Phys. D: Appl. Phys. 33 1817 (2000).
  * [12] “Generation of a cylindrically symmetric, polarized laser beam with narrow linewidth and fine tunability,” T. Hirayama, Y. Kozawa, T. Nakamura, and S. Sato, Opt. Express 14, 12839 (2006).
  * [13] “Radially and azimuthally polarized beams generated by space-variant dielectric subwavelength gratings,” Z. Bomzon, G. Biener, V. Kleiner, and E. Hasman, Opt. Lett. 27, 285 (2002).
  * [14] “Generation of a radially polarized laser beam by use of the birefringence of a c-cut Nd:YVO4 crystal,” K. Yonezawa, Y. Kozawa, and S. Sato, Opt. Lett. 31, 2151 (2006).
  * [15] “A versatile and stable device allowing the efficient generation of beams with radial, azimuthal or hybrid polarizations,” T. Grosjean, A. Sabac and D. Courjon, Opt. Commun. 252, 12 (2005)
  * [16] “Generation of the rotationally symmetric TE01and TM01modes from a wavelength-tunable laser,” J. J. Wynne, IEEE J. Quantum Electron. 10, 125 (1974).
  * [17] “Generation of a radially polarized laser beam by use of a conical Brewster prism,” Y. Kozawa and S. Sato, Opt. Lett. 30, 3063 (2005).
  * [18] “Fields of a radially polarized Gaussian laser beam beyond the paraxial approximation,” Y. I. Salamin, Opt. Lett. 31, 2619 (2006).
  * [19] “2kW, M2 $<$ 10 radially polarized beams from aberration-compensated rod-based Nd:YAG lasers,” I. Moshe, S. Jackel, A. Meir, Y. Lumer, and E. Leibush, Opt. Lett. 32, 47 (2007).
  * [20] “Production of radially or azimuthally polarized beams in solid-state lasers and the elimination of thermally induced birefringence effects,” I. Moshe, S. Jackel, and A. Meir, Opt. Lett. 28, 807 (2003).
  * [21] “Generation of radially polarized beams in a Nd:YAG laser with self-adaptive overcompensation of the thermal lens,” M. Roth, E. Wyss, H. Glur, and H.P. Weber Opt. Lett. 30, 1665 (2005).
  * [22] “Influence of beam polarization on laser cutting efficency,” V. G. Niziev, A. V. Nestorov, J. Phys. D: Appl. Phys. 32, 1455, (1999).
  * [23] “Linear, annular, and radial focusing with axicons and applications to laser machining,” M. Rioux, R. Tremblay, and P.-A. Belanger, Appl. Opt. 17, 1532 (1978).
  * [24] “Mono-energetic GeV electrons from ionization in a radially polarized laser beam,” Y. I. Salamin, Opt. Lett. 32, 90 (2007).
  * [25] “A tunable doughnut laser beam for cold-atom experiments,” S. M. Iftiquar, J. Opt. B: Quantum Semiclass. Opt. 5, 40 (2003).
  * [26] “Focusing of high numerical aperture cylindrical-vector beams,” K. Youngworth and T. Brown, Opt. Express 7, 77 (2000).
  * [27] “Sharper focus for a radially polarized beam,” R. Dorn, S. Quabis, and G. Leuches Phys. Rev. Lett. 91, 233901 (2003).
  * [28] “An all-fiber device for generating radially and other polarized light beams,” T. Grosjean, D. Courjon and M. Spajer, Opt. Commun. 203, 1 (2002).
  * [29] “Generation of cylindrical vector beams with few-mode fibers excited by Laguerre Gaussian beams,” G. Volpe and D. Petrov, Opt. Commun. 237, 89 (2004).
  * [30] “Generation of radially polarized mode in Yb fiber laser by using a dual conical prism,” J. Li, K. Ueda, M. Musha, A. Shirakawa, and L. Zhong, Opt. Lett. 31, 2969 (2006).


#  Dihedral Lattice Gauge Theories on a Quantum Annealer

Michael Fromm<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  Christopher Winterowd
<EMAIL_ADDRESS>Institut für Theoretische Physik, Goethe-
Universität Frankfurt
Max-von-Laue-Str. 1, 60438 Frankfurt am Main, Germany

###### Abstract

We study lattice gauge theory with discrete, non-Abelian gauge groups. We
extend the formalism of previous studies on D-Wave’s quantum annealer as a
computing platform to finite, simply reducible gauge groups. As an example, we
use the dihedral group $D_{n}$ with $n=3,4$ on a two plaquette ladder for
which we provide proof-of-principle calculations of the ground-state and
employ the known time evolution formalism with Feynman clock states.

Hamiltonian lattice gauge theory, Quantum computing, Quantum annealing

##  I Introduction

Lattice gauge theory (LGT) calculations using quantum computers have already
seen substantial progress. This is despite the fact that programmable quantum
hardware has only recently become widely available to researchers in physics
(see e.g. [1] for an up-to-date high-energy physics perspective on the field).
With the basic formulation of quantum simulations of LGT being laid out very
early on [2], efficient formulations of abelian [3] and non-abelian LGT [4, 5]
on universal, gate-based hardware now exist. This includes a complete set of
instructions for the efficient and accurate simulation of QCD and QED [6].

On the side of adiabatic quantum computing [7], despite the fact that it has
been commercially available for more than a decade in the form of quantum
annealers (QA) [8], this approach has only recently been used for LGT
calculations. These include the pioneering studies on the annealer for the
case of $SU(2)$ [9] and $SU(3)$ [10]. In these formulations, the number of
qubits necessary to digitize the theory under study scales with the size of
the Hilbert space of the problem rather than the spatial volume of the system.
Thus, this formulation does not show the expected quantum advantage present in
universal quantum computing. On the other hand, systems in a QA architecture
s.a. D-Wave’s Advantage_system5.1 already comprise several thousand physical
qubits [11]. Even at this stage of hardware development, proof-of-principle
quantum computations in LGT [9, 10] or other field theories [12] are feasible.
The intrinsic nature of the formulation of problems on the annealer simply
requires the mapping of the lattice field theory onto an optimization problem
represented by the Hamiltonian

$\displaystyle H(q)=\sum_{i}Q_{ii}q_{i}+\sum_{i<j}Q_{ij}q_{i}q_{j}\,,$ (1)

with a real, upper-triangular matrix $Q$ and binary variables
$q_{i}\in\\{0,1\\}$. This type of problem goes under the name of quadratic
unconstrained binary optimization (QUBO). Thus, QA can be seen as an ideal
entry point into the field of quantum simulation of LGT. It is in this spirit,
that we turn to the subject of our study.

On the lattice, the gauge field Hilbert space is infinite for compact Lie
groups and a truncation of the full symmetry becomes necessary (see [13] for a
basic overview). One such truncation is the approximation of the full symmetry
group by a discrete subgroup, a strategy whose varied success and intricacies
are well summarized in [14]. Here, we take the practical approach by choosing
the finite, non-abelian dihedral group $D_{n}$ as our gauge group, which can
be digitized without truncation as in [15]. By mapping this problem onto an
optimatization problem amenable to the QA, we extend the previous studies of
compact Lie groups [9, 10] to the case of simply reducible, finite groups for
which we provide the adapted framework.

##  II Hamiltonian $D_{n}$ lattice gauge theory

We begin by introducing the Hamiltonian formulation of $D_{n}$ lattice gauge
theory. Our approach closely follows that of [2], where the Hamiltonian
approach was worked out for a general gauge group $G$. As usual, we work on a
cubic lattice of dimension $d$, with the gauge fields living on the links
between the lattice sites. The gauge field Hilbert space, denoted by
$\mathcal{H_{G}}$, is a direct product of the individual link spaces
$\mathcal{H_{G}}=\bigotimes_{\ell}\mathcal{H}_{\ell}$ which, in the group
element basis, are defined by
$\mathcal{H}_{\ell}=\mathrm{span}(\\{\ket{g}\\}_{g\in G})$. We start by
defining the link operator $\hat{U}$, acting on $\mathcal{H}_{\ell}$, as

$\displaystyle\hat{U}^{j}_{mn}=\int\mathrm{dg}\,D^{j}_{mn}(g)\ket{g}\bra{g}$
(2)

where $j$ labels the irreducible representations (irreps) of $G$ and
$D^{j}_{mn}(g)$ are the Wigner representation matrices for $g\in G$. The
indices $m,n$ label the multiplicity of the left and right projection of the
link, respectively. This object is of primary importance in the Hamiltonian
formulation as it is responsible for the interactions.

Figure 1:  Ladder geometry used in this study. The sites of the ladder have
been labeled as $(i,j)$, where $i=0,1,\dots,N-1$ and $j=0,1$. Here $N$ denotes
the length of the ladder. The forward link operators are also shown and have
been labeled by their direction $\mu=1,2$ as well as the site from which they
emanate. We note that the system is periodic only in the $\hat{1}$-direction.

The Hamiltonian formulation for $G=SU(2)$ lattice gauge theory was first
written down in [16]. In later work, it was shown how this Hamiltonian could
be obtained from the transfer matrix [17]. For a general gauge group $G$, the
lattice Hamiltonian, commonly referred to as the Kogut-Susskind Hamiltonian,
consists of two terms

$\displaystyle\hat{H}_{\text{KS}}=\hat{H}_{\text{E}}+\hat{H}_{\text{B}},$ (3)

where the first term is referred to as the electric part and the second term
is referred to as the magnetic part. The magnetic part takes the form

$\displaystyle\hat{H}_{\text{B}}=-\lambda_{B}\sum_{\vec{x},i<j}\mathrm{Re}\operatorname{Tr}\hat{U}_{i}(\vec{x})\hat{U}_{j}(\vec{x}+\hat{i})\hat{U}^{\dagger}_{i}(\vec{x}+\hat{j})\hat{U}^{\dagger}_{j}(\vec{x}),$

where the sum is over the spatially-oriented plaquettes111For the case of
$G=D_{3},D_{4}$ which we will consider in this work, we take the link
operators to be in the faithful, two-dimensional irrep, which we denote by
$\mathbb{2}$. and can be seen to be diagonal in the group elements basis by
virtue of the definition in Eq.(2). In this basis, the electric term,
$\hat{H}_{\text{E}}$, has a form that depends on whether $G$ is taken to be a
compact Lie group or a finite group [18, 19]. It is often much more convenient
to work in the representation basis, labeled by the states $\ket{jmn}$. Here
$j$ labels the irrep and $m,n$ run over the states within the multiplet. In
this basis $\hat{H}_{\text{E}}$ becomes diagonal

$\displaystyle\hat{H}_{\text{E}}=\lambda_{E}\sum_{\vec{x}}\sum^{d}_{i=1}\sum_{jmn}f_{j}\ket{jmn}_{\vec{x},i}\bra{jmn}_{\vec{x},i}\,.$
(5)

The group element basis is related to the representation basis through the
following relation

$\displaystyle\braket{g}{jmn}=\sqrt{\frac{\text{dim}(j)}{|G|}}D^{j}_{mn}(g),$
(6)

where $|G|$ is the order of the finite group $G$. From here on out, we will
assume that we are working with finite groups and our formulae will reflect
this. Using the transformation in Eq.(6), one can easily transform between the
two bases.

It is important to discuss the couplings in the Kogut-Susskind Hamiltonian.
The relationship between the couplings $\lambda_{E},\lambda_{B}$ appearing in
front of the electric and magnetic terms, respectively, can be determined in
the limit $a_{0}\to 0$ when deriving the Hamiltonian in the transfer matrix
formalism. However, this procedure depends on the gauge group under
consideration. While for a compact Lie group,
$\lambda_{E}=g_{H}^{2}/2,\lambda_{B}=1/g_{H}^{2}$, for a finite group it has
been determined by previous studies that one should use
$\lambda_{E}=\exp{(-2/g^{2}_{H})},\lambda_{B}=1/g_{H}^{2}$, instead [18, 19].
Here we have introduced the Hamiltonian coupling $g_{H}$, which is the
geometric mean of the spatial and temporal couplings, $g_{s}$ and $g_{t}$.
These are introduced in the transfer matrix formulation when one takes the
temporal and spatial lattice spacings to be distinct. The coefficients $f_{j}$
appearing in Eq.(5) are the eigenvalues of the quadratic Casimir operator for
the case of compact Lie groups [2]. For example, when $G=SU(2)$, one gets the
familiar result $f_{j}=j(j+1)$. On the other hand, for finite gauge groups,
the $f_{j}$ can be derived systematically from the transfer matrix in the
limit of vanishing temporal lattice spacing, $a_{0}\to 0$ [19]. Other choices
for the coefficients for the case of discrete gauge groups also exist in the
literature [20].

###  II.1 Computation of $H_{ij}$

With the Kogut-Susskind Hamiltonian given in operator form Eq.(3), there still
remains the task of constructing the states of the physical Hilbert space
$\mathcal{H_{P}}$. The label physical here refers to the subspace of the
larger Hilbert space $\mathcal{H_{G}}$ that respects local gauge invariance.
Formally, in the group element basis, local gauge invariance can be expressed
with the help of left and right multiplication operators given by
$\Theta^{L}_{g}(\vec{x},i)$ and $\Theta^{R}_{g}(\vec{x},i)$ that act as
follows

$\displaystyle\Theta^{L}_{g}\ket{h}=\ket{g^{-1}h},\leavevmode\nobreak\
\Theta^{R}_{g}\ket{h}=\ket{hg^{-1}},\leavevmode\nobreak\ g,h\in G,$ (7)

with $\Theta^{R\dagger}_{g}=\Theta^{R}_{g^{-1}}$. As usual, the link
transforms in the adjoint representation under a gauge transformation. Thus, a
local gauge transformation at the site $\vec{x}$ parametrized by $g$ is given
by

$\displaystyle\tilde{\Theta}_{g}(\vec{x})\equiv\prod^{d}_{i=1}\Theta^{L}_{g}(\vec{x},i)\Theta^{R\dagger}_{g}(\vec{x}-\hat{i},i).$
(8)

For a generic physical state $\ket{\psi}$, gauge invariance demands

$\displaystyle\tilde{\Theta}_{g}(\vec{x})\ket{\psi}=\ket{\psi},\forall\leavevmode\nobreak\
\vec{x}\in\mathbb{Z}^{d}.$ (9)

In the representation basis, the statement of Gauss’s law in Eq.(9) is
equivalent to $\ket{\psi}$ being written as a direct product of color singlets
at each lattice site

$\displaystyle\ket{\psi}=\bigotimes_{\vec{x}}\ket{00}_{\vec{x}}\,,$ (10)

where we refer to Appendix B for the details rearding the explicit
construction of $\ket{00}_{\vec{x}}$ in terms of the $\ket{jmn}$.

As we are ultimately interested in mapping our system onto a quantum annealer,
we are restricted to a rather small system size (Fig. (1)). For a given
lattice geometry, we are then left with the task of determining the matrix
elements of the Hamiltonian (3). We start by introducing the trivial vacuum
state $\ket{0}$ with every link $\ell$ in the trivial representation,
$j_{\ell}=0$, s.t. $\hat{H}_{\text{E}}\ket{0}=0$. Physically, this corresponds
to the situation where the links contain zero chromoelectric flux. Following
the approach of [2], the states of the physical Hilbert space can be
systematically generated from this configuration by subsequently acting with
gauge-invariant operators. To carry out this procedure, one needs to know how
an individual link operator acts on a general state in the representation
basis, $\ket{jmn}$. Using the definition of the link operator in (2) and the
matrix element in (6), one obtains

$\displaystyle\hat{U}^{\mathbb{2}}_{m^{\prime}n^{\prime}}\ket{jmn}=\sqrt{\frac{\text{dim}(j)}{|G|}}\sum_{g\in
G}D^{2}_{m^{\prime}n^{\prime}}(g)D^{j}_{mn}(g)\ket{g}\,.$

This can be further simplified by using the Clebsch-Gordan (CG) series for the
tensor product of two arbitrary representation matrices, which yields the
result

$\displaystyle\hat{U}^{\mathbb{2}}_{m^{\prime}n^{\prime}}\ket{jmn}$
$\displaystyle=$
$\displaystyle\sum_{J}\sum_{M}\sum_{N}\sqrt{\frac{\text{dim}(j)}{\text{dim}(J)}}$
$\displaystyle\times$
$\displaystyle\braket{2m^{\prime}jm}{JM}\braket{JM}{2n^{\prime}jn}\ket{JMN},$

where the sum on $J$ is over the irreps and the sums over $M$ and $N$ are over
the states in a given irrep.222We note that for our choice of the $\mathbb{2}$
irrep of $D_{3}$ and $D_{4}$ Eq.(25), the link operators satisfy
$U^{\mathbb{2}\dagger}_{m,n}=U^{\mathbb{2}}_{n,m}$.

Applying gauge invariant operators repeatedly to the trivial vacuum state with
the help of Eq.(II.1), the enumeration of the configuration space can be
performed. In this procedure, Gauss’s law is imposed at each step. This is
equivalent to the Wigner $3$J-symbol being nonzero for a tuple of irreps,
$(j_{1},j_{2},j_{3})$, which characterize the three links involving a given
lattice site. For more details regarding the $3$J symbols and Gauss’s law we
refer the reader to appendices (A.3) and (B).

This task of mapping out the physical Hilbert space can be automated using a
Markov-chain-like approach. For this, all one needs are the $3$J symbols for
the given gauge group $G$. The total number of states in the full Hilbert
space is $\tilde{N}_{\text{irreps}}^{3N}$ states, where
$\tilde{N}_{\text{irreps}}$ is the total number of irreps of the gauge group
$G$ and $N$ is the size of the ladder. Although enforcing local gauge
invariance removes a large number of these states, the physical Hilbert space
still grows quite rapidly. In Table 1, we display the size of the physical
Hilbert space, where $G=D_{3},D_{4}$. These counts show that even for modest
system sizes and small non-abelian groups there are constraints as to the
problems that can be mapped to the quantum annealer.

$N$ | 2 | 3 | 4
---|---|---|---
$D_{3}$ | 49 | 251 | $O(1300)$
$D_{4}$ | 76 | 392 | $O(2500)$
Table 1: List of the size of the physical Hilbert space, $N_{\text{conf}}$, on
a ladder of size $N$ for $D_{3}$ and $D_{4}$. The configurations are
enumerated by a set of integers $\\{j_{i},i=1,2,\dots,3N\\}$ characterizing
the irrep of each link on whereby Gauss’s law is satisfied at each site.

Once we have enumerated the states in the physical Hilbert, space we can
finally compute the matrix elements of the Hamiltonian in this basis. As
$\hat{H}_{\text{E}}$ is diagonal in the representation basis, the application
of Eq.(5) to the generic state in (84) is trivial. The magnetic Hamiltonian,
which is responsible for the interactions, has a nontrivial action on a
physical state. To illustrate this, we first write an arbitrary plaquette
operator on the ladder

$\displaystyle\hat{P}_{\vec{x}}$ $\displaystyle=$
$\displaystyle\sum_{n_{1},\ldots,n_{4}}\hat{U}_{\vec{x},1;n_{1},n_{2}}\hat{U}_{\vec{x}+\hat{1},2;n_{2},n_{3}}$
(13) $\displaystyle\times$
$\displaystyle\hat{U}_{\vec{x}+\hat{2},1;n_{3},n_{4}}^{\dagger}\hat{U}_{\vec{x},1;n_{4},n_{1}}^{\dagger},$

where $\vec{x}=(x,0),\leavevmode\nobreak\ x=0,\ldots,N-1$ denotes the vertex
at the bottom left corner of the plaquette and the sum is over the group
indices. Here we label each link by its site vector and direction. Now, using
the relation Eq.(II.1), we can determine the result of the plaquette operator
acting on a state. The matrix element of Eq.(13) between two arbitrary
physical states is given by

$\displaystyle\bra{\psi^{\prime}}\hat{P}_{\vec{x}}\ket{\psi}$ $\displaystyle=$
$\displaystyle\sum_{n_{1},\ldots,n_{4}}\sum_{\\{m_{i,1},o_{i,1}\\}}\cdots\sum_{\\{m_{i,n_{s}},o_{i,n_{s}}\\}}$
$\displaystyle\times$
$\displaystyle\prod_{s}\begin{pmatrix}j_{\ell_{1,s}}&j_{\ell_{2,s}}&j_{\ell_{3,s}}\\\
m_{1,s}&m_{2,s}&m_{3,s}\end{pmatrix}\overline{\begin{pmatrix}l_{\ell_{1,s}}&l_{\ell_{2,s}}&l_{\ell_{3,s}}\\\
o_{1,s}&o_{2,s}&o_{3,s}\end{pmatrix}}$ $\displaystyle\times$
$\displaystyle\prod_{i\not\in\mathcal{L}_{\vec{x}}}\delta_{l_{i},j_{i}}\delta_{m_{R_{i}},o_{R_{i}}}\delta_{m_{L_{i}},o_{L_{i}}}$
$\displaystyle\times$
$\displaystyle\frac{\sqrt{\mathrm{dim}(j_{\tilde{l}_{1}})\mathrm{dim}(j_{\tilde{l}_{2}})\mathrm{dim}(j_{\tilde{l}_{3}})\mathrm{dim}(j_{\tilde{l}_{4}})}}{\sqrt{\mathrm{dim}(l_{\tilde{l}_{1}})\mathrm{dim}(l_{\tilde{l}_{2}})\mathrm{dim}(l_{\tilde{l}_{3}})\mathrm{dim}(l_{\tilde{l}_{4}})}}$
$\displaystyle\times$
$\displaystyle\braket{2n_{1}j_{\tilde{l}_{1}}m_{L_{\tilde{l}_{1}}}}{l_{\tilde{l}_{1}}o_{L_{\tilde{l}_{1}}}}\braket{l_{\tilde{l}_{1}}o_{R_{\tilde{l}_{1}}}}{2n_{2}j_{\tilde{l}_{1}}m_{R_{\tilde{l}_{1}}}}$
$\displaystyle\times$
$\displaystyle\braket{2n_{2}j_{\tilde{l}_{2}}m_{L_{\tilde{l}_{2}}}}{l_{\tilde{l}_{2}}o_{L_{\tilde{l}_{2}}}}\braket{l_{\tilde{l}_{2}}o_{R_{\tilde{l}_{2}}}}{2n_{3}j_{\tilde{l}_{2}}m_{R_{\tilde{l}_{2}}}}$
$\displaystyle\times$
$\displaystyle\braket{2n_{4}j_{\tilde{l}_{3}}m_{L_{\tilde{l}_{3}}}}{l_{\tilde{l}_{3}}o_{L_{\tilde{l}_{3}}}}\braket{l_{\tilde{l}_{3}}o_{R_{\tilde{l}_{3}}}}{2n_{3}j_{\tilde{l}_{3}}m_{R_{\tilde{l}_{3}}}}$
$\displaystyle\times$
$\displaystyle\braket{2n_{1}j_{\tilde{l}_{4}}m_{L_{\tilde{l}_{4}}}}{l_{\tilde{l}_{4}}o_{L_{\tilde{l}_{4}}}}\braket{l_{\tilde{l}_{4}}o_{R_{\tilde{l}_{4}}}}{2n_{4}j_{\tilde{l}_{4}}m_{R_{\tilde{l}_{4}}}},$

where we refer to $\ket{\psi}$ and $\ket{\psi^{\prime}}$ as the “in" and “out"
states,
$\mathcal{L}_{\vec{x}}=\\{\tilde{l}_{1},\tilde{l}_{2},\tilde{l}_{3},\tilde{l}_{4}\\}$
denotes the set of all links involved in (13), and the bar denotes complex
conjugation. The sums over $n_{i}$ run over the states in the $\mathbb{2}$
representation and the sums over $m_{i,\mu},o_{i,\mu}$ run over the states in
the corresponding irreps for each link in the “in" and “out" states. In
Eq.(II.1), we have introduced the Wigner $3$J symbols which can be generalized
to $D_{n}$. This along with other details regarding the group theory of
$D_{n}$ are given in appendix (A), with the full derivation of Eq.(II.1) in
appendix (C). Using this result for the plaquette matrix element, one can
construct the magnetic Hamiltonian by recalling that
$\hat{H}_{\text{B}}=-\lambda_{B}\sum_{\vec{x}}\hat{P}_{\vec{x}}$. We note here
that this construction is completely general and applies to a ladder of
arbitrary length.

Figure 2:  Visualization of the Hamiltonian matrices ($g_{H}^{2}=0.75$) for
gauge groups $D_{3}$ (left) and $D_{4}$ (right) on the $N=2$ ladder. One
immediately notices the sparsity of both Hamiltonians, which is due to the
structure of the magnetic contribution.

For the gauge groups which we have examined, it turns out that the Hamiltonian
matrix is extremely sparse. This will work to our advantage later on when we
map our problem to the quantum annealer. In Fig.2, we display a visualization
of the sparsity of the Hamiltonian for both $D_{3}$ and $D_{4}$. Examining the
product of delta functions in Eq.(II.1), the sparseness of the Hamiltonian
ultimately comes from the orthonormality of the link states in the
representation basis. Once one has calculated the Hamiltonian matrix, the
classical part of the calculation is practically complete. In the following,
we discuss how our lattice gauge Hamiltonian is transformed into an
optimization problem so that the quantum computation on the annealer can be
performed.

##  III Implementation and Results

###  III.1 Groundstate via variational formulation

Quantum annealing is a method used to solve a very specific type of problem:
the calculation of the ground state of a generalized Ising model [11]. Thus,
unlike the case of the gate-based approach where one has at one’s disposal a
set of universal quantum gates, for quantum annealing one must cast the
problem that one would like to solve into the form of an Ising model.

In the context of lattice gauge theory, it has been shown that one can map the
Hamiltonian in Eq.(3) onto a model with QUBO form, Eq.(1), in order to compute
the low-lying states of the spectrum [10, 9]. To see how this emerges we
consider the variational principle from quantum mechanics

$\displaystyle
E_{0}\leq\frac{\bra{\psi}\hat{H}\ket{\psi}}{\braket{\psi}{\psi}}\,,$ (15)

where $E_{0}$ is the ground-state energy. Here we use a variational ansatz
with a trial wave function

$\displaystyle\ket{\psi}=\sum_{\alpha=1}^{N_{\mathrm{conf}}}a_{\alpha}\ket{\phi_{\alpha}},$
(16)

with real parameters $a_{\alpha}$ and basis states $\ket{\phi_{\alpha}}$
corresponding to the configurations of the physical Hilbert space. The
expansion parameters are, in general, complex but here can be chosen to be
real as the Hamiltonian is a real, symmetric matrix. It is in this way that
solving for the ground-state energy of our lattice Hamiltonian can be recast
as an optimization problem as one seeks to minimize
$\bra{\psi}\hat{H}\ket{\psi}$. Our accuracy in determining the eigenstates of
the system is limited only by the precision to which we can determine the
coefficients $a_{\alpha}$. To cast our problem in the form of Eq.(1), we do
the following. First, the coefficients are given a fixed-point binary
representation. Second, the norm of the wave function is discouraged from
being zero by adding a penalty term. Incorporating both of these into our
variational calculation, the rhs of Eq.(15) can rewritten as a cost function

$\displaystyle
F=\bra{\psi}\hat{H}\ket{\psi}-\eta\braket{\psi}{\psi}=\sum_{\alpha,\beta}^{{N_{\mathrm{conf}}}}\sum_{i,j}^{K}Q_{\alpha\beta,ij}q_{\alpha,i}q_{\beta,j},$

Figure 3: (_Left_): Ground-state expectation values for $D_{3}$ for the full
Kogut-Susskind Hamiltonian (red), electric term (black), and magnetic term
(blue) as a function of the inverse Hamiltonian coupling squared. The open
symbols represent the minimum result from the quantum annealer while the
filled symbols were obtained from classical simulated annealing. The colored
band, although barely visible, represents the mean and sample standard
deviation of the measurements from the quantum annealer. Simulation parameters
for the latter were $K=3$ with $z_{max}=5$ zoom steps and
$n_{\mathrm{reads}}=1000$. (_Right_): The same quantities for $G=D_{4}$ where
the sample standard deviation from QA is much larger. This is due to the fact
that more computing resources are needed to accurately determine the minimum.
Simulation parameters for QA were $K=2$ with $z_{max}=7$ zoom steps and
$n_{\mathrm{reads}}=2000$.

where

$\displaystyle Q_{\alpha\beta,ij}$ $\displaystyle=$ $\displaystyle
2^{i+j-2K-2z}(-1)^{\delta_{iK}+\delta_{jK}}h_{\alpha\beta}+\delta_{\alpha\beta}\delta_{ij}\tilde{Q}_{\alpha,i}$
$\displaystyle\tilde{Q}_{\alpha,i}$ $\displaystyle=$ $\displaystyle
2^{i-K-z+1}(-1)^{\delta_{iK}}\sum_{\gamma}^{N_{conf}}a_{\gamma}^{(z)}h_{\gamma\alpha}\,$
(18) $\displaystyle h_{\alpha\beta}$ $\displaystyle=$ $\displaystyle
H_{\alpha\beta}-\eta\delta_{\alpha\beta}.$

Here, $\eta$ represents the tunable parameter multiplying the penalty term. In
addition to giving the variational parameters a floating-point representation,
in the above definitions we have already introduced the parameter $z$, which
is used in the adaptive variational search method [10]. This procedure
iteratively improves the estimates for the $a^{(z+1)}_{\alpha}$ by
distributing the $K$ sampling points around the preceding solution to the Qubo
problem, $a^{(z)}_{\alpha}$,

$\displaystyle a_{\alpha}^{(z+1)}$ $\displaystyle=$ $\displaystyle
a_{\alpha}^{(z)}-\frac{q_{\alpha,K}}{2^{z}}+\sum_{i=1}^{K-1}\frac{q_{\alpha,i}}{2^{K-i+z}}\,,$
(19)

starting at $a_{\alpha}^{(0)}=0$. The estimate for the eigenstate is refined
at each step, hence the name “zooming" for this procedure. The number of zoom
steps plays a significant role in the overall computational cost of our
calculation, as each refinement requires calls to the quantum annealer.

The quantum computations are done on the quantum annealing hardware
Advantage_system5.1 from D-Wave [11], which is accessible via its API D-Wave
Ocean [21]. As our system sizes are still small (c.f. Table 1), the results
can be compared with the exact solution using the Hamiltonian Eq.(3) as well
as with simulated annealing via the Ocean package neal. As alternative, which
we did not employ, D-Wave offers hybrid solvers s.a. the KerberosSampler [21]
that attempt to break down the original QUBO matrix into smaller pieces to be
subsequently solved using classical or quantum hardware. This appears to be
particularly useful for system sizes that cannot be embedded on currently
available annealer architectures.

It should be noted that for computations employing both, simulated annealing
and quantum annealing, results still have a $\eta$ dependence, see Eq.(III.1).
As already noticed in [10, 9], convergence to the true ground-state is
achieved for $\eta$ lying in the actual vicinity of the ground-state energy,
$E_{0}$ (approaching from above). For practical reasons, one can determine the
“suitable" $\eta$ for a given $Q$ by solving Eq.(III.1) iteratively in $\eta$,
terminating the calculation when a certain convergence criterion, s.a.
relative improvement in the solution, is fulfilled. This strategy works well
for local computations with simulated annealing and could in principle be
employed also for quantum annealing. Here, runtime on the quantum annealer is
the major constraint.

We finally comment on our setup when accessing the annealer via the provided
software package. We use the quantum annealer in its forward annealing mode
with default annealing schedule and annealing time
$t_{f}=20\leavevmode\nobreak\ \mathrm{\mu s}$. At least one more parameter
needs to be provided by the user during the quantum annealing computations.
This is what is referred to as the chain strength. For our calculations, we
find automatic chain strength tuning (default option) to be sufficient. Fig.3
shows results from both, simulated annealing and quantum annealing, for the
ground-state energy $\langle H\rangle$ (red) as well as the expectation value
for the magnet part $\langle H_{B}\rangle$ (blue) and kinetic part $\langle
H_{E}\rangle$ (black) for $G=D_{3}$ (left). When going to $G=D_{4}$ (right),
an increase in computational resources is needed due to the more complicated
energy landscape for the larger group, which we however only partially meet
due to runtime restrictions.

###  III.2 Time Evolution

One of the main motivations for working in the Hamiltonian formulation of
lattice gauge theories is the ability to access real-time dynamics. This
stands in stark contrast to mainstream lattice calculations which work in
Euclidean space and must perform an analytic continuation of numerical data in
order to access real-time physics. In the gate-based approach, the so-called
Trotter approximation is employed to the time evolution operator,
$\hat{U}(T)=\exp\\{-iT\hat{H}\\}$, which evolves an initial state by a finite
time $T$ [15]. This approximates consists of replacing the full $\hat{U}$ with
products of operators which evolve the system on a smaller time interval,
$\delta t$. Corrections to this approximation typically scale with powers of
$\delta t$. This approach to time-evolution of quantum states allows for an
efficient simulation of the theory using universal quantum computers.

In order to solve this problem on the quantum annealer, however, one must
reformulate time evolution as an optimization problem. This can be done by the
introduction of Feynman clock states [22], a mechanism first applied to
quantum chemistry calculations in order to generate parallel-in-time quantum
dynamics [23]. We thus have to introduce an ancillary quantum system with
states $\ket{t},t=1,2,\dots,N_{t}$ where $N_{t}$ is the number of time-slices
in the time evolution. Tensoring this orthonormal state with our as-yet-
unknown state vector $\ket{\psi_{t}}$ at each timeslice, the problem of time
evolution is equivalent to the minimization of the following functional

$\displaystyle\mathcal{L}$ $\displaystyle=$
$\displaystyle\sum^{N_{t}}_{t,t^{\prime}=1}\bra{t^{\prime}}\bra{\psi_{t^{\prime}}}\hat{\mathcal{C}}\ket{\psi_{t}}\ket{t}$
(20) $\displaystyle-$
$\displaystyle\eta\left(\sum^{N_{t}}_{t,t^{\prime}=1}\bra{t^{\prime}}\braket{\psi_{t^{\prime}}}{\psi_{t}}\ket{t}-1\right),$

where

$\displaystyle\hat{\mathcal{C}}$ $\displaystyle\equiv$
$\displaystyle\hat{\mathcal{C}}_{0}+\frac{1}{2}\bigg{(}\mathbb{I}\otimes\ket{t}\bra{t}+\mathbb{I}\otimes\ket{t+1}\bra{t+1}$
(21) $\displaystyle-$ $\displaystyle\hat{U}_{\delta
t}\otimes\ket{t+1}\bra{t}-\hat{U}_{\delta
t}^{\dagger}\otimes\ket{t}\bra{t+1}\bigg{)}.$

Figure 4:  Results for time evolution using simulated annealing at
$g_{H}^{2}=0.75$ with simulation parameters given in the text. The red and
yellow points represent the expectation value of the magnetic Hamiltonian in
the time-evolved trivial vacuum state as a function of time for $D_{3}$ and
$D_{4}$. The blue and black points represent the probability amplitude for the
trivial vacuum state to persist as a function of time. The lines represent the
exact result for each case.

Here $\delta t\equiv T/N_{t}$ is the step size in time, $\eta$ is a Lagrange
multiplier analogous to our previous penalty term, and $\hat{\mathcal{C}}_{0}$
selects a predetermined initial state. By construction, $\hat{\mathcal{C}}$ is
hermitian. One can show that the minimum of the functional Eq.(20) corresponds
to the exact time-evolved state at each step. Thus, as a result of a single
optimization problem one obtains the the full time-evolution of a many-body
quantum state over a finite time interval. From this functional one can now
obtain the QUBO matrix. As previously discussed for the case of finding the
ground state of our Hamiltonian, this is what the quantum annealer requires as
input. Our discussion closely follows the derivation of [24]. Using a
variational state $\ket{\psi_{t}}\ket{t}$ at each time step $t$, the
functional in Eq.(20) becomes

$\displaystyle\mathcal{L}=\sum_{\alpha\beta}a^{*}_{\alpha}L_{\alpha\beta}a_{\beta},$
(22)

where the expansion parameters $a_{\alpha}$ are complex, the indices in the
sum run over all $N_{t}N_{\mathrm{conf}}$ values, and $L_{\alpha\beta}$ are
the matrix elements of the functional. The terms in the above sum can be
written in terms of the real and imaginary parts of both $L_{\alpha\beta}$ and
$a_{\alpha}$. Using the fixed-point representation for both the real and
imaginary parts of $a_{\alpha}$, one obtains the following QUBO matrix

$Q_{\alpha,i;\beta,j}=\begin{cases}2^{i+j-2K-2z}(-1)^{\delta_{iK}+\delta_{jK}}\Re
L_{\alpha\beta}+2\delta_{\alpha\beta}\delta_{ij}2^{i-K-z}(-1)^{\delta_{iK}}\sum_{\gamma}\left(\Re
a^{(z)}_{\gamma}\Re L_{\gamma\beta}+\Im a^{(z)}_{\gamma}\Im
L_{\gamma\beta}\right)&1\leq i,j\leq K\,,\\\
-2^{i+j^{\prime}-2K-2z}(-1)^{\delta_{iK}+\delta_{j^{\prime}K}}\Im
L_{\alpha\beta}&1\leq i,j^{\prime}\leq K\,,\\\
2^{i^{\prime}+j-2K-2z}(-1)^{\delta_{i^{\prime}K}+\delta_{jK}}\Im
L_{\alpha\beta}&1\leq i^{\prime},j\leq K,,\\\
2^{i^{\prime}+j^{\prime}-2K-2z}(-1)^{\delta_{i^{\prime}K}+\delta_{j^{\prime}K}}\Re
L_{\alpha\beta}+2\delta_{\alpha\beta}\delta_{i^{\prime}j^{\prime}}2^{i^{\prime}-K-z}(-1)^{\delta_{i^{\prime}K}}\sum_{\gamma}\left(\Im
a^{(z)}_{\gamma}\Re L_{\gamma\beta}-\Re a^{(z)}_{\gamma}\Im
L_{\gamma\beta}\right)&1\leq i^{\prime},j^{\prime}\leq K\,.\end{cases}$ (23)

where the Latin indices now run from $1$ to $2K$ to allow for $K$ bits in
representing both the real and imaginary parts of the variational parameters
and the primed Latin indices are shifted by $K$. The dimension of this QUBO
matrix is $2KN_{t}N_{\mathrm{conf}}\times 2KN_{t}N_{\mathrm{conf}}$, which is
significantly larger than the one used to determine the eigenstates of the
Hamiltonian. This problem size is too large to be fully embedded on current
quantum annealers and we revert to local simulated annealing333Simulation
parameters are $K=3,N_{t}=3,z_{\mathrm{max}}=11$ and $n_{\mathrm{runs}}=8000$.
to solve it. The results of our time evolution simulations are displayed in
Fig.(4) where we have time evolved the trivial vacuum. Shown are the
expectation value of the magnetic part, $\hat{H}_{B}$, as well as the
probability that the trivial vacuum persists,
$\left|\bra{0}U(t)\ket{0}\right|^{2}$. One can see good agreement with the
exact results.

##  IV Conclusion and Outlook

We have demonstrated how to map a simple, non-abelian lattice gauge theory
onto a quantum annealer. In doing so, we have been able to compute the
spectrum as well as time evolution for both $D_{3}$ and $D_{4}$ on small
lattices. Furthermore, we have provided a detailed explanation of how to
construct the physical states for the Kogut-Susskind Hamiltonian for an
arbitrary group. In principle, this construction could be generalized to both
larger groups as well as larger lattices. The main obstacle in doing so is the
scaling of the physical Hilbert space with increasing lattice and group size.
It thus may be helpful to look in other directions in order to utilize the
power of the quantum annealer.

One idea to get around the barrier of rapidly expanding Hilbert spaces for the
case of continuous groups is to perform a truncation in the number of allowed
irreps as was done in earlier studies of $SU(2)$ and $SU(3)$ in the
Hamiltonian formulation [4, 5, 9]. This could also be done for $D_{n}$, where
$n>4$. A similar truncation procedure would involve both the electric and
magnetic terms in the Hamiltonian and thus one could investigate the affects
on the ground-state energy as well as the dynamics of the system.

One further avenue that could be pursued with the annealer is the estimation
of tunneling rates and vacuum decay [25, 12]. These non-perturbative processes
are fundamental to the understanding of a wide range of phenomena in both
high-energy and condensed matter physics. Another possibility is using quantum
annealing in state preparation. This is an important problem faced by gate-
based approaches to simulating lattice gauge theories.

###### Acknowledgements.

The authors would like to thank Brendan Pankovich for useful discussions and
Michael Spannowsky for comments on the manuscript. 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).

## References

  * [1] 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, “Quantum simulation for high energy physics,” Apr. 2022.
  * [2] T. Byrnes and Y. Yamamoto, “Simulating lattice gauge theories on a quantum computer,” Phys.Rev. A73 (2006) 022328, Oct. 2005.
  * [3] J. F. Haase, L. Dellantonio, A. Celi, D. Paulson, A. Kan, K. Jansen, and C. A. Muschik, “A resource efficient approach for quantum and classical simulations of gauge theories in particle physics,” Quantum 5, 393 (2021), June 2020.
  * [4] N. Klco, J. R. Stryker, and M. J. Savage, “Su(2) non-abelian gauge field theory in one dimension on digital quantum computers,” Phys. Rev. D 101, 074512 (2020), Aug. 2019.
  * [5] A. Ciavarella, N. Klco, and M. J. Savage, “A trailhead for quantum simulation of su(3) yang-mills lattice gauge theory in the local multiplet basis,” Phys. Rev. D 103, 094501 (2021), Jan. 2021.
  * [6] A. Kan and Y. Nam, “Lattice quantum chromodynamics and electrodynamics on a universal quantum computer,” July 2021.
  * [7] T. Albash and D. A. Lidar, “Adiabatic quantum computation,” Rev. Mod. Phys., vol. 90, p. 015002, Jan 2018.
  * [8] M. W. Johnson, M. H. S. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk, E. M. Chapple, C. Enderud, J. P. Hilton, K. Karimi, E. Ladizinsky, N. Ladizinsky, T. Oh, I. Perminov, C. Rich, M. C. Thom, E. Tolkacheva, C. J. S. Truncik, S. Uchaikin, J. Wang, B. Wilson, and G. Rose, “Quantum annealing with manufactured spins,” Nature, vol. 473, pp. 194–198, may 2011.
  * [9] S. A. Rahman, R. Lewis, E. Mendicelli, and S. Powell, “Su(2) lattice gauge theory on a quantum annealer,” Phys. Rev. D 104, 034501 (2021), Mar. 2021\.
  * [10] M. Illa and M. J. Savage, “Basic elements for simulations of standard model physics with quantum annealers: Multigrid and clock states,” Feb. 2022.
  * [11] https://www.dwavesys.com, D-Wave Systems Inc.
  * [12] S. Abel and M. Spannowsky, “Observing the fate of the false vacuum with a quantum laboratory,” PRX Quantum 2, 010349 (2021), Mar. 2021.
  * [13] E. Zohar, “Quantum simulation of lattice gauge theories in more than one space dimension – requirements, challenges, methods,” June 2021.
  * [14] M. S. Alam, S. Hadfield, H. Lamm, and A. C. Y. Li, “Quantum simulation of dihedral gauge theories,” Aug. 2021.
  * [15] H. Lamm, S. Lawrence, and Y. Yamauchi, “General methods for digital quantum simulation of gauge theories,” Phys. Rev. D 100, 034518 (2019), Mar. 2019\.
  * [16] J. Kogut and L. Susskind, “Hamiltonian formulation of wilson’s lattice gauge theories,” Phys. Rev. D, vol. 11, pp. 395–408, Jan 1975.
  * [17] M. Creutz, “Gauge fixing, the transfer matrix, and confinement on a lattice,” Phys. Rev. D, vol. 15, pp. 1128–1136, Feb 1977.
  * [18] D. Harlow and H. Ooguri, “Symmetries in quantum field theory and quantum gravity,” 2018.
  * [19] M. Carena, H. Lamm, Y.-Y. Li, and W. Liu, “Lattice renormalization of quantum simulations,” Phys. Rev. D, vol. 104, p. 094519, Nov 2021.
  * [20] J. Bender, E. Zohar, A. Farace, and J. I. Cirac, “Digital quantum simulation of lattice gauge theories in three spatial dimensions,” New J. Phys. 20, 093001 (2018), Apr. 2018.
  * [21] https://docs.ocean.dwavesys.com.
  * [22] R. P. Feynman, “Quantum mechanical computers,” Optics News, vol. 11, pp. 11–20, Feb 1985.
  * [23] J. R. McClean, J. A. Parkhill, and A. Aspuru-Guzik, “Feynman’s clock, a new variational principle, and parallel-in-time quantum dynamics,” Proceedings of the National Academy of Sciences, vol. 110, no. 41, pp. E3901–E3909, 2013.
  * [24] M. Illa and M. J. Savage, “Basic elements for simulations of standard model physics with quantum annealers: Multigrid and clock states,” 2022.
  * [25] S. Abel, N. Chancellor, and M. Spannowsky, “Quantum computing for quantum tunneling,” Phys. Rev. D, vol. 103, p. 016008, Jan 2021.
  * [26] E. P. Wigner, A. S. Wightman, and J. Mehra, The collected works of Eugene Paul Wigner.  Spinger-Verlag, 1992-1998.
  * [27] P. M. van Den Broek and J. F. Cornwell, “Clebsch-gordan coefficients of symmetry groups,” physica status solidi (b), vol. 90, no. 1, pp. 211–224, 1978.
  * [28] Y.-T. Hu, Emergent properties in exactly solvable discrete models for two-dimensional topological phases.  PhD thesis, University of Utah, 2013.
  * [29] Y. Hu, N. Geer, and Y.-S. Wu, “Full dyon excitation spectrum in extended levin-wen models,” Phys. Rev. B, vol. 97, p. 195154, May 2018.

##  Appendix A The Gauge Group $D_{n}$

In this appendix we summarize key facts and properties of the dihedral group,
$D_{n}$, that are relevant to this work. In particular, we list the
irreducible representations for general $n$, the Clebsch-Gordan coefficients
for both $D_{3}$ and $D_{4}$, and outline the algorithm for constructing the
Wigner $3$J-symbols.

###  A.1 Irreducible Representations

One can define the dihedral group as the symmetry group of a regular,
$n$-sided polygon. The symmetry operations are thus, rotation and reflections
in the plane in which the polygon resides. The order of the group, denoted by
$|D_{n}|$, is $2n$. We denote the group element which corresponds to rotations
by $2\pi/n$ by $r$ and that which corresponds to a reflection about, say, the
$y$-axis by $s$. These two elements satisfy the following relations

$\displaystyle r^{n}=s^{2}=e,\leavevmode\nobreak\ srs=r^{-1}=r^{n-1},$ (24)

where $e$ denotes the identity element. Any element of $D_{n}$ can be uniquely
represented either as $r^{k},0\leq k<n$, or as $sr^{k},0\leq k<n$. We now list
the irreducible representations (irreps) of $D_{n}$ for arbitrary $n$. We
recall that a representation of a group is irreducible if it contains no
invariant subspaces.

We start with the one-dimensional irreps of the dihedral group. For $n$ even,
there are four one-dimensional irreps while for $n$ odd, there are two one-
dimensional irreps. This is summarized by the characters listed in Table 2.
Two-dimensional irreps of $D_{n}$ for arbitrary $n$ also exist and are given,
in the standard basis, by

$\displaystyle\rho^{h}_{s}(r^{j})=\begin{pmatrix}\cos\frac{2\pi
hj}{n}&-\sin\frac{2\pi hj}{n}\\\ \sin\frac{2\pi hj}{n}&\cos\frac{2\pi
hj}{n}\end{pmatrix},\ \ \rho^{h}_{s}(sr^{j})=\begin{pmatrix}-\cos\frac{2\pi
hj}{n}&\sin\frac{2\pi hj}{n}\\\ \sin\frac{2\pi hj}{n}&\cos\frac{2\pi
hj}{n}\end{pmatrix},$ (25)

where we require $0<h<n/2$, where $h$ labels the various irreps. From (25),
the characters can be easily read off

$\displaystyle\chi_{h}(r^{j})=2\cos\left(\frac{2\pi
hj}{n}\right),\leavevmode\nobreak\ \chi_{h}(sr^{j})=0.$ (26)

These will prove useful when decomposing an arbitrary tensor representation
into a direct sum of irreps.

| $r^{k}$ | $sr^{k}$
---|---|---
$\chi_{1}(g)$ | 1 | 1
$\chi_{2}(g)$ | 1 | -1
$\chi_{3}(g)$ | $(-1)^{k}$ | $(-1)^{k}$
$\chi_{4}(g)$ | $(-1)^{k}$ | $(-1)^{k+1}$
Table 2: Character table for the one-dimensional irreps of $D_{n}$. For $n$
even, all four irreps are valid, while for $n$ odd, only the first two rows
apply. In discussing $D_{3}$, we will commonly refer to these one-dimensional
irreps as the trivial representation, denoted by $\mathbb{1}$, and the
$\sigma$-representation.

###  A.2 Clebsch-Gordan Coefficients

An important task in many branches of physics is to understand how one
decomposes a tensor product of two irreps into a direct sum of irreps. For the
case of $D_{n}$ lattice gauge theory, in order to construct the physical
Hilbert space by acting with gauge-invariant loop operators, one must perform
exactly this task. We consider two unitary irreps of a group $G$, $\rho^{p}$
and $\rho^{q}$. In general, the tensor product of these two representations,
denoted by $\rho^{p}\otimes\rho^{q}$, is reducible. This implies that it is
equivalent to a direct sum of unitary irreps, which can be symbolically
expressed as

$\displaystyle\rho^{p}\otimes\rho^{q}\approx\sum_{r}\oplus\leavevmode\nobreak\
n^{r}_{p,q}\rho^{r},$ (27)

where $n^{r}_{p,q}$ is the number of times that the irrep labeled by $r$
appears in the direct sum. It is immediately clear that $n^{r}_{p,q}$ is
symmetric in $p$ and $q$. The dihedral group $D_{n}$ is an example of a so-
called simply reducible group, where, among other things, $n_{r}=0,1$ [26].
The following analysis is made simpler by the requirement of only working with
simply reducible groups, although it can be generalized to groups where the
multiplicities can take values larger than one [27].

The Clebsch-Gordan (CG) coefficients can be thought of as a non-singular
matrix transformation $U$, from the direct product basis to the block-diagonal
basis. This is expressed by the following relation

$\displaystyle\left(\rho^{p}\otimes\rho^{q}\right)U=U\left(\sum_{r}\oplus\leavevmode\nobreak\
n^{r}_{p,q}\rho^{r}\right),$ (28)

there $U$ is of the size $d_{p}d_{q}\times d_{p}d_{q}$, where $d_{r}$ denotes
the size of the irrep $r$. The CG coefficients are denoted by

$\displaystyle U_{i,j;r,\alpha,k}\equiv\left(\begin{array}[]{cc|c}p&q&r\\\
i&j&k\end{array}\right),$ (31)

where $i=1,2,\dots,d_{p}$, $j=1,2,\dots,d_{q}$, and $k=1,2,\dots,d_{r}$. We
list here the results for $D_{3}$ and $D_{4}$, where in principle, the CG
coefficients can be determined for arbitrary $n$.

####  A.2.1 $D_{3}$

The tensor products of the one-dimensional irreps are straightforward to
compute as
$n^{\sigma}_{\mathbb{1},\sigma}=n^{\mathbb{1}}_{\sigma,\sigma}=n^{\mathbb{1}}_{\mathbb{1},\mathbb{1}}=1$.
We list below the results for the tensor products involving the two-
dimensional irrep (labeled as $\mathbb{2}$)

$\displaystyle\left(\begin{array}[]{cc|c}2&\mathbb{1}&2\\\
m&&l\end{array}\right)=\delta_{ml},$ (34)
$\displaystyle\left(\begin{array}[]{cc|c}2&\sigma&2\\\
m&&l\end{array}\right)=\epsilon_{ml},$ (37)
$\displaystyle\left(\begin{array}[]{cc|c}2&2&\mathbb{1}\\\
m&n&\end{array}\right)=\frac{\delta_{mn}}{\sqrt{2}},$ (40)
$\displaystyle\left(\begin{array}[]{cc|c}2&2&\sigma\\\
m&n&\end{array}\right)=\frac{\epsilon_{mn}}{\sqrt{2}},$ (43)
$\displaystyle\left(\begin{array}[]{cc|c}2&2&2\\\
m&n&l\end{array}\right)=\frac{1}{\sqrt{2}}\left(\delta_{l2}\sigma^{z}_{mn}-\delta_{l1}\sigma^{x}_{mn}\right),$
(46)

where the indices $m,n,l=1,2$, enumerating the basis elements of the two-
dimensional irrep.

####  A.2.2 $D_{4}$

With the addition of two additional one-dimensional irreps, the case of
$D_{4}$ is a straightforward extension of the calculations for $D_{3}$. The
tensor products of the one-dimensional irreps are summarized by

$\displaystyle
n^{2}_{34}=n^{4}_{23}=n^{3}_{24}=n^{q}_{1q}=1,\leavevmode\nobreak\
\tilde{q}=1,2,3,4,$ (47)

where we use the same labeling of the one-dimensional irreps as in Table 2.
The CG coefficients for the $\mathbb{2}$ (two-dimensional) irrep tensored with
the various one-dimensional irreps are as follows

$\displaystyle\left(\begin{array}[]{cc|c}\mathbb{2}&2&\mathbb{2}\\\
m&&l\end{array}\right)=\epsilon_{ml},$ (50)
$\displaystyle\left(\begin{array}[]{cc|c}\mathbb{2}&3&\mathbb{2}\\\
m&&l\end{array}\right)=\sigma^{z}_{ml},$ (53)
$\displaystyle\left(\begin{array}[]{cc|c}\mathbb{2}&4&\mathbb{2}\\\
m&&l\end{array}\right)=\sigma^{x}_{ml},$ (56)
$\displaystyle\left(\begin{array}[]{cc|c}\mathbb{2}&1&\mathbb{2}\\\
m&&l\end{array}\right)=\delta_{ml}.$ (59)

The tensor product of the $\mathbb{2}$ irrep with itself gives a direct sum of
all four one-dimensional irreps. The CG coefficients associated with this
decomposition are given by

$\displaystyle\left(\begin{array}[]{cc|c}\mathbb{2}&\mathbb{2}&1\\\
m&n&\end{array}\right)=\frac{\delta_{mn}}{\sqrt{2}},$ (62)
$\displaystyle\left(\begin{array}[]{cc|c}\mathbb{2}&\mathbb{2}&2\\\
m&n&\end{array}\right)=\frac{\epsilon_{mn}}{\sqrt{2}},$ (65)
$\displaystyle\left(\begin{array}[]{cc|c}\mathbb{2}&\mathbb{2}&3\\\
m&n&\end{array}\right)=\frac{\sigma^{z}_{mn}}{\sqrt{2}},$ (68)
$\displaystyle\left(\begin{array}[]{cc|c}\mathbb{2}&\mathbb{2}&4\\\
m&n&\end{array}\right)=\frac{\sigma^{x}_{mn}}{\sqrt{2}}.$ (71)

For $D_{n},n>4$, the analysis becomes more involved as multiple two-
dimensional irreps appear.

###  A.3 Wigner $3$J-symbols

In order to construct the physical Hilbert space we need to obtain the Wigner
$3$J-symbols for $D_{n}$. Although the construction for $SU(2)$ is familiar to
most, it turns out that there exists a generalization to finite groups [28,
29].

Before providing the relevant formulae, we first provide a few key definitions
from representation theory which will aid our discussion. We start by
considering a generic irrep of a finite group denoted by $(\rho,V)$, where
$\rho$ is a homomorphism from the group elements to a vector space $V$. To
each irrep $j$ we can associate a dual rep, denoted by $j^{*}$, which is
equivalent to the rep $\\{\rho^{*}(g),\leavevmode\nobreak\ \forall g\in G\\}$.
The irrep $j$ is related to its dual via the so-called duality map
$\Omega^{j}$, which is unitary and commutes with the action of the
representation matrices. Furthermore, the duality map has the property that
$(\Omega^{j})^{T}=\alpha_{j}\Omega^{j}$, where $\alpha_{j}=\pm 1,0$ is the
Frobenius-Schur (FS) indicator.

We denote the 3$J$ symbol for the irreps of group $G$ labeled by the tuple
$(j_{1},j_{2},j_{3})$ by

$\displaystyle
C_{j_{1}j_{2}j_{3}m_{1}m_{2}m_{3}}\equiv\begin{pmatrix}j_{1}&j_{2}&j_{3}\\\
m_{1}&m_{2}&m_{3}\end{pmatrix},$ (72)

where the right-hand side uses the familiar notation from $SU(2)$. The cyclic
property of the $3$J symbol is given by

$\displaystyle
C_{j_{1}j_{2}j_{3};m_{1}m_{2}m_{3}}=\alpha_{j_{3}}C_{j_{3}j_{1}j_{2};m_{3}m_{1}m_{2}},$
(73)

while one can obtain the $3$J symbol for the dual of the tuple of irreps using

$\displaystyle
C_{j^{*}_{3}j^{*}_{2}j^{*}_{1};m_{3}m_{2}m_{1}}=\sum_{n_{1}n_{2}n_{3}}\overline{C_{j_{1}j_{2}j_{3};n_{1}n_{2}n_{3}}\Omega^{j^{*}_{3}}_{n_{3}m_{3}}\Omega^{j^{*}_{2}}_{n_{2}m_{2}}\Omega^{j^{*}_{1}}_{n_{1}m_{1}}},$
(74)

where the bar represents complex conjugation. Finally, the normalization of
the $3$J symbol is fixed by requiring

$\displaystyle\sum_{\\{m_{i}\\},\\{n_{i}\\}}\Omega^{j_{3}}_{m_{3}n_{3}}\Omega^{j_{2}}_{m_{2}n_{2}}\Omega^{j_{1}}_{m_{1}n_{1}}C_{j_{1}j_{2}j_{3};m_{1}m_{2}m_{3}}C_{j^{*}_{3}j^{*}_{2}j^{*}_{1};n_{3}n_{2}n_{1}}=1.$
(75)

This guarantees that when we couple three irreps together at a site of our
ladder via the $3$J symbol, the resulting singlet state will be normalized to
one. The duality map as well as the $3$J symbol can be constructed for
arbitrary irreps by using the algorithm presented in [28].

###  A.4 $D_{3}$

The results for the $3$J symbols for $D_{3}$ are as follows

$\displaystyle
C_{\mathbb{1}\mathbb{1}\mathbb{1}}=C_{\mathbb{1}\sigma\sigma}=1,$ (76)
$\displaystyle
C_{\mathbb{1}\mathbb{2}\mathbb{2};m_{1}m_{2}}=\begin{pmatrix}\frac{1}{\sqrt{2}}&0\\\
0&\frac{1}{\sqrt{2}}\end{pmatrix}_{m1,m2},$ (77) $\displaystyle
C_{\sigma\mathbb{2}\mathbb{2};m_{1}m_{2}}=\begin{pmatrix}0&\frac{i}{\sqrt{2}}\\\
\frac{-i}{\sqrt{2}}&0\end{pmatrix}_{m1,m2},$ (78) $\displaystyle
C_{\mathbb{2}\mathbb{2}\mathbb{2};m_{1}m_{2}m_{3}}=\begin{pmatrix}\\{0,\frac{1}{2}\\}&\\{\frac{1}{2},0\\}\\\
\\{\frac{1}{2},0\\}&\\{0,-\frac{1}{2}\\}\end{pmatrix}_{m1,m2,m_{3}}.$ (79)

###  A.5 $D_{4}$

In a similar fashion, one can obtain the $3$J symbols for $D_{4}$, which read

$\displaystyle
C_{111}=C_{122}=C_{133}=C_{144}=C_{1\mathbb{2}\mathbb{2}}=C_{34\mathbb{2}}=1,$
(80) $\displaystyle
C_{1\mathbb{2}\mathbb{2}}=\begin{pmatrix}\frac{1}{\sqrt{2}}&0\\\
0&\frac{1}{\sqrt{2}}\end{pmatrix}_{m1,m2},C_{2\mathbb{2}\mathbb{2}}=\begin{pmatrix}0&\frac{i}{\sqrt{2}}\\\
-\frac{i}{\sqrt{2}}&0\end{pmatrix}_{m1,m2},$ (81) $\displaystyle
C_{3\mathbb{2}\mathbb{2}}=\begin{pmatrix}\frac{1}{\sqrt{2}}&0\\\
0&-\frac{1}{\sqrt{2}}\end{pmatrix}_{m1,m2},C_{2\mathbb{2}\mathbb{2}}=\begin{pmatrix}0&\frac{1}{\sqrt{2}}\\\
\frac{1}{\sqrt{2}}&0\end{pmatrix}_{m1,m2}.$ (82)

##  Appendix B Gauss’s Law

For the construction of the physical Hilbert space, one must enforce Gauss’s
law at each lattice site. This is precisely the expression of local gauge
invariance. For the more common case where compact Lie groups are employed,
this condition can be written down in terms of the generators of left and
right gauge transformations on the links which enter or leave a given site of
the lattice. An analogous construction, which can be generalized to discrete
groups such as $D_{n}$, uses properties of the states characterizing the
various irreps. To be more concrete, at each site of our ladder we tensor
together the three links which either start or end on that site while summing
over the internal degrees of freedom such that it is a singlet. This has a
natural expression in terms of the Wigner $3$J-symbols, and can be written as
follows

$\displaystyle\ket{00}$ $\displaystyle=$
$\displaystyle\sum_{m_{\text{I}},\tilde{m}_{\text{I}}}\ket{j_{\text{I}},m_{\text{I}};j_{\text{I}},\tilde{m}_{\text{I}}}\braket{j_{\text{I}},m_{\text{I}};j_{\text{I}},\tilde{m}_{\text{I}}}{0,0},$
(83) $\displaystyle=$
$\displaystyle\sum_{m_{\text{I}},\tilde{m}_{\text{I}}}\sum_{m_{\text{A}}}\sum_{m_{\text{E}}}\ket{j_{\text{A}}m_{\text{A}}}\otimes\ket{j_{\text{E}}m_{\text{E}}}\otimes\ket{j_{\text{I}}m_{\text{I}}}\braket{j_{\text{A}},m_{\text{A}};j_{\text{E}},m_{\text{E}}}{j_{\text{I}},\tilde{m}_{\text{I}}}\braket{j_{\text{I}},m_{\text{I}};j_{\text{I}},\tilde{m}_{\text{I}}}{0,0},$
$\displaystyle=$
$\displaystyle\sum_{m_{\text{I}},\tilde{m}_{\text{I}}}\sum_{m_{\text{A}}}\sum_{m_{\text{E}}}(-)^{f(j_{\text{A}},j_{\text{E}},j_{\text{I}};\tilde{m}_{\text{I}})}\ket{j_{\text{A}}m_{\text{A}}}\otimes\ket{j_{\text{E}}m_{\text{E}}}\otimes\ket{j_{\text{I}}m_{\text{I}}}\begin{pmatrix}j_{\text{A}}&j_{\text{E}}&j_{\text{I}}\\\
m_{\text{A}}&m_{\text{E}}&\tilde{m}_{\text{I}}\end{pmatrix}\braket{j_{\text{I}},m_{\text{I}};j_{\text{I}},\tilde{m}_{\text{I}}}{0,0},$
$\displaystyle=$
$\displaystyle\sum_{m_{\text{I}}}\sum_{m_{\text{A}}}\sum_{m_{\text{E}}}(-)^{f(j_{\text{A}},j_{\text{E}},j_{\text{I}};m_{\text{I}})}\ket{j_{\text{A}}m_{\text{A}}}\otimes\ket{j_{\text{E}}m_{\text{E}}}\otimes\ket{j_{\text{I}}m_{\text{I}}}\begin{pmatrix}j_{\text{A}}&j_{\text{E}}&j_{\text{I}}\\\
m_{\text{A}}&m_{\text{E}}&m_{\text{I}}\end{pmatrix},$

where $f$ is a known function which has been tabulated for $D_{3}$ and $D_{4}$
and takes integer or half-integer values. In the third line we have used the
fact that the CG coefficients for the projection of the tensor product of any
two irreps onto the trivial representation ($\mathbb{1}$) is diagonal. The
expression in (83) is valid for continuous as well as finite groups. Using the
Wigner $3$J-symbols and the CG coefficients for $D_{3}$ and $D_{4}$ which were
listed in appendix (A), one can produce a gauge invariant state on the ladder
by taking the tensor product of the analogous expression at each site. This is
represented by the expression

$\displaystyle\ket{\psi}=\bigotimes_{s}\ket{00}_{s}=\sum_{\\{m_{i,1}\\}}\ldots\sum_{\\{m_{i,n_{s}}\\}}\left(\bigotimes_{\ell}\ket{j_{\ell}m_{L_{\ell}}m_{R_{\ell}}}\right)\prod_{s}\begin{pmatrix}j_{\ell_{1,s}}&j_{\ell_{2,s}}&j_{\ell_{3,s}}\\\
m_{1,s}&m_{2,s}&m_{3,s}\end{pmatrix}(-)^{f(j_{\ell_{1,s}},j_{\ell_{2,s}},j_{\ell_{3,s}};m_{3,s})},$
(84)

where the $\\{m_{i,s}\\}$ are the local multiplicities at site $s$ of the link
$\ell$ in representation $j_{\ell}$, and in the product over links, each link
has multiplicities $m_{L_{\ell}}$ and $m_{R_{\ell}}$ associated with the
“left" and “right" ends of the link. One can thus label the states of the
Hilbert space by a set of $N_{\text{links}}=3N$ integers characterizing the
irreps, where $N$ is the number of plaquettes, with the condition that at each
site there is a valid tuple with nonzero $3$J symbol. From here, one has
everything one needs in order to enumerate the states in the physical Hilbert
space.

##  Appendix C Plaquette Contribution

Here we provide the details for the matrix element of a plaquette in the
physical Hilbert space. This derivation is similar to that of [9] for the case
of $SU(2)$. The application of the plaquette in (13) on the arbitrary physical
state listed in (84) is given by

$\displaystyle\hat{P}_{\vec{x}}\ket{\psi}$ $\displaystyle=$
$\displaystyle\sum_{n_{1},\ldots,n_{4}}\sum_{\\{m_{L_{i}}\\}}\sum_{\\{m_{R_{i}}\\}}\hat{U}_{\tilde{l}_{1};n_{1}n_{2}}\ket{j_{\tilde{l}_{1}}m_{L_{\tilde{l}_{1}}}m_{R_{\tilde{l}_{1}}}}$
$\displaystyle\otimes$
$\displaystyle\hat{U}_{\tilde{l}_{2};n_{2}n_{3}}\ket{j_{\tilde{l}_{2}}m_{L_{\tilde{l}_{2}}}m_{R_{\tilde{l}_{2}}}}\otimes\hat{U}_{\tilde{l}_{3};n_{4}n_{3}}\ket{j_{\tilde{l}_{3}}m_{L_{\tilde{l}_{3}}}m_{R_{\tilde{l}_{3}}}}$
$\displaystyle\otimes$
$\displaystyle\hat{U}_{\tilde{l}_{4};n_{1}n_{4}}\ket{j_{\tilde{l}_{4}}m_{L_{\tilde{l}_{4}}}m_{R_{\tilde{l}_{4}}}}\bigotimes_{i\not\in\mathcal{L}_{\vec{x}}}\ket{j_{i}m_{L_{i}}m_{R_{i}}}$
$\displaystyle\prod_{s}$
$\displaystyle\begin{pmatrix}j_{\ell_{1,s}}&j_{\ell_{2,s}}&j_{\ell_{3,s}}\\\
m_{1,s}&m_{2,s}&m_{3,s}\end{pmatrix}(-)^{f(j_{\ell_{1,s}},j_{\ell_{2,s}},j_{\ell_{3,s}};m_{3,s})}\,,$
(86)

where the sums over $m_{L_{i}}$ and $m_{R_{i}}$ represent $3N$ individual sums
whose ranges depend on the configuration of the state $\ket{\psi}$. Using the
Clebsch-Gordan series, we can obtain the result of each link operator acting
on an arbitrary link state. Applying this to each term in (86) yields

$\displaystyle\hat{P}_{\vec{x}}\ket{\psi}=\sum_{n_{1},\ldots,n_{4}}\sum_{\\{m_{L_{i}}\\}}\sum_{\\{m_{R_{i}}\\}}\prod_{s}\begin{pmatrix}j_{\ell_{1,s}}&j_{\ell_{2,s}}&j_{\ell_{3,s}}\\\
m_{1,s}&m_{2,s}&m_{3,s}\end{pmatrix}(-)^{f(j_{\ell_{1,s}},j_{\ell_{2,s}},j_{\ell_{3,s}};m_{3,s})}\bigotimes_{i\not\in\mathcal{L}_{\vec{x}}}\ket{j_{i}m_{L_{i}}m_{R_{i}}}\otimes$
$\displaystyle\sum_{J_{\tilde{l}_{1}}}\sum_{M_{L_{\tilde{l}_{1}}}}\sum_{M_{R_{\tilde{l}_{1}}}}\sqrt{\frac{\mathrm{dim}(j_{\tilde{l}_{1}})}{\mathrm{dim}(J_{\tilde{l}_{1}})}}\braket{2n_{1}j_{\tilde{l}_{1}}m_{L_{\tilde{l}_{1}}}}{J_{\tilde{l}_{1}}M_{L_{\tilde{l}_{1}}}}\braket{J_{\tilde{l}_{1}}M_{R_{\tilde{l}_{1}}}}{2n_{2}j_{\tilde{l}_{1}}m_{R_{\tilde{l}_{1}}}}\ket{J_{\tilde{l}_{1}}M_{L_{\tilde{l}_{1}}}M_{R_{\tilde{l}_{1}}}}\otimes$
$\displaystyle\sum_{J_{\tilde{l}_{2}}}\sum_{M_{L_{\tilde{l}_{2}}}}\sum_{M_{R_{\tilde{l}_{2}}}}\sqrt{\frac{\mathrm{dim}(j_{\tilde{l}_{2}})}{\mathrm{dim}(J_{\tilde{l}_{2}})}}\braket{2n_{2}j_{\tilde{l}_{2}}m_{L_{\tilde{l}_{2}}}}{J_{\tilde{l}_{2}}M_{L_{\tilde{l}_{2}}}}\braket{J_{\tilde{l}_{2}}M_{R_{\tilde{l}_{2}}}}{2n_{3}j_{\tilde{l}_{2}}m_{R_{\tilde{l}_{2}}}}\ket{J_{\tilde{l}_{2}}M_{L_{\tilde{l}_{2}}}M_{R_{\tilde{l}_{2}}}}\otimes$
$\displaystyle\sum_{J_{\tilde{l}_{3}}}\sum_{M_{L_{\tilde{l}_{3}}}}\sum_{M_{R_{\tilde{l}_{3}}}}\sqrt{\frac{\mathrm{dim}(j_{\tilde{l}_{3}})}{\mathrm{dim}(J_{\tilde{l}_{3}})}}\braket{2n_{4}j_{\tilde{l}_{3}}m_{L_{\tilde{l}_{3}}}}{J_{\tilde{l}_{3}}M_{L_{\tilde{l}_{3}}}}\braket{J_{\tilde{l}_{3}}M_{R_{\tilde{l}_{3}}}}{2n_{3}j_{\tilde{l}_{3}}m_{R_{\tilde{l}_{3}}}}\ket{J_{\tilde{l}_{3}}M_{L_{\tilde{l}_{3}}}M_{R_{\tilde{l}_{3}}}}\otimes$
$\displaystyle\sum_{J_{\tilde{l}_{4}}}\sum_{M_{L_{\tilde{l}_{4}}}}\sum_{M_{R_{\tilde{l}_{4}}}}\sqrt{\frac{\mathrm{dim}(j_{\tilde{l}_{4}})}{\mathrm{dim}(J_{\tilde{l}_{4}})}}\braket{2n_{1}j_{\tilde{l}_{4}}m_{L_{\tilde{l}_{4}}}}{J_{\tilde{l}_{4}}M_{L_{\tilde{l}_{4}}}}\braket{J_{\tilde{l}_{4}}M_{R_{\tilde{l}_{4}}}}{2n_{4}j_{\tilde{l}_{4}}m_{R_{\tilde{l}_{4}}}}\ket{J_{\tilde{l}_{4}}M_{L_{\tilde{l}_{4}}}M_{R_{\tilde{l}_{4}}}}\,.$
(87)

We now take the inner product of (87) with a second generic state
$\ket{\psi^{\prime}}$ which satisfies Gauss’ law. We note that states in the
representation basis are orthonormal

$\displaystyle\braket{j^{\prime}m^{\prime}n^{\prime}}{jmn}=\delta_{jj^{\prime}}\delta_{mm^{\prime}}\delta_{nn^{\prime}},$
(88)

and thus we obtain the result in Eq.(II.1).