#  Overcoming exponential scaling with system size in Trotter-Suzuki
implementations of constrained Hamiltonians: 2+1 U(1) lattice gauge theories

Dorota M. Grabowska<EMAIL_ADDRESS>Theoretical Physics Department,
CERN, 1211 Geneva 23, Switzerland  Christopher Kane<EMAIL_ADDRESS>Department of Physics, University of Arizona, Tucson, AZ 85719, USA  Benjamin
Nachman<EMAIL_ADDRESS>Physics Division, Lawrence Berkeley National
Laboratory, Berkeley, CA 94720, USA  Christian W. Bauer<EMAIL_ADDRESS>Physics Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720,
USA

###### Abstract

For many quantum systems of interest, the classical computational cost of
simulating their time evolution scales exponentially in the system size. At
the same time, quantum computers have been shown to allow for simulations of
some of these systems using resources that scale polynomially with the system
size. Given the potential for using quantum computers for simulations that are
not feasible using classical devices, it is paramount that one studies the
scaling of quantum algorithms carefully. This work identifies a term in the
Hamiltonian of a class of constrained systems that naively requires quantum
resources that scale exponentially in the system size. An important example is
a compact ${\rm U}(1)$ gauge theory on lattices with periodic boundary
conditions. Imposing the magnetic Gauss’ law a priori introduces a constraint
into that Hamiltonian that naively results in an exponentially deep circuit. A
method is then developed that reduces this scaling to polynomial in the system
size, using a redefinition of the operator basis. An explicit construction of
the matrices defining the change of operator basis, as well as the scaling of
the associated computational cost, is given.

††preprint: CERN-TH-2022-133

##  I Introduction

Constrained Hamiltonians are ubiquitous in fundamental physics, as any gauge
theory – general relativity, electromagnetism, Yang–Mills, string theory – all
result in such Hamiltonians. Unfortunately, constrained Hamiltonians pose a
number of difficulties, both theoretical and practical. On the theoretical
side, it is not always clear how to quantize the theory such that all
constraints are properly implemented in a self-consistent manner [1, 2, 3]. On
the practical side, as will be discussed in this paper, constrained
Hamiltonians can result in a large degree of coupling, resulting in resource-
expensive simulations.

This work will mostly be concerned with the “degree of coupling” (DoC) between
various terms in the Hamiltonian. The DoC of a given term is defined as the
maximum number of operators that need to be included in a single term of a
Suzuki-Trotter represenation [4, 5, 6]. Furthermore, the DoC of a Hamiltonian
is defined as the maximum DoC of any given term within the Hamiltonian itself.
For example, for a Hamiltonian of the form
$H=\sum^{N_{\mathcal{O}}}_{ij}\mathcal{O}_{i}\mathcal{O}_{j}$, where the
operators $\mathcal{O}_{i}$ are local operators and $N_{\mathcal{O}}$ is the
total number of such operators, the ${\rm DoC}=2$, independent of the total
number of operators. On the other hand, a Hamiltonian of the form
$H=F(\sum_{i}^{N_{\mathcal{O}}}\mathcal{O}_{i})$ will have

$\displaystyle{\rm
DoC}[F;N_{\mathcal{O}}]\equiv\min\left[N_{\mathcal{O}},\deg(F)\right]$ (1)

where $\deg(F)$ is the polynomial degree of the function $F$; note that a non-
polynomial function has $\deg(F)=\infty$. Throughout this work, the arguments
of a given DoC will sometimes be suppressed for brevity.

As we will discuss, constraints often give rise to Hamiltonians where the DoC
scales with the system size, which in turn naively creates an exponential
scaling in the gates required to simulate the time evolution of such a system.
Note that the degree of coupling is different from the non-locality of the
system, which typically has to do with the minimum distance between two
entangled qubits. A well-known, and phenomenologically relevant, class of
constrained Hamiltonians are gauge theories formulated on space-time lattices
with periodic boundary conditions. In such theories, the Hilbert space is
divided into different charge sectors, with the gauge invariant Hamiltonian
not including any interactions that migrate between the gauge sectors.
However, due to noisy machines as well as a practical inability to exactly
simulate time evolution, one needs to contend with the different charge
sectors. Formulating these theories using only the physical subspace
introduces the constraints with which this paper is concerned. It is important
to note, however, that there are proposals for quantum simulations that do not
rely on constructing Hamiltonians that span only the physical Hilbert space.
The basic idea of these methods is to introduce an energy penalty term in the
Hamiltonian that suppresses transitions that are gauge variant [7, 8, 9]. For
more details on this approach see Refs. [10, 11, 12, 13, 7, 14, 15, 16, 17,
18].

As a test case, the dual-basis formulation of a ${\rm U}(1)$ lattice gauge
theory in two spatial dimensions111For the remainder of this manuscript, we
will use ‘D+1’ dimensions to refer to a theory that has D spatial dimensions
and one time dimension. is studied [19, 20, 21, 22]. While this theory does
not appear in nature, it is a very useful case study. In particular, the
compact formulation is expected to display linear confinement, similar to QCD
in three spatial dimensions [23, 24, 25, 26, 27, 19]. As discussed in Ref.
[20], for such a theory, the set of magnetic fields on each plaquette of the
lattice is generally over-complete, since the magnetic Gauss’ law constrains
the total magnetic flux through any closed surface to vanish. For a 2+1
dimensional lattice with periodic boundary conditions, this leads to one
constraint. For open boundary conditions, there is no closed surface and so no
constraint appears. In 3+1 dimensions, the number of constraints for periodic
boundary conditions is $2+N$, where $N$ is defined as the total number of
lattice sites; in this dimension, it is also the total number of cubes in the
volume [20]. The number of constraints can be understood by realizing that
there are $N-1$ independent constraints on single-cell cubes and then three
extensive constraints involving a quadratic number of plaquette operators, one
for each plaquette orientation (in the $\hat{x}-\hat{y}$ plane,
$\hat{x}-\hat{z}$ plane and $\hat{y}-\hat{z}$ plane); note that these three
extensive constraints do not appear for open boundary conditions. Solving the
constraints and writing one of the operators in a given closed surface in
terms of the other operators in that same surface can remedy the redundancy of
the basis. Doing so, however, introduces terms in the Hamiltonian that couple
together many operators at different lattice sites, naively giving rise to
exponentially scaling circuits on a quantum device.

This manuscript presents a general method for breaking this exponential
scaling in system size, based on the idea of reducing the degree of coupling
via a judiciously chosen change in the operator basis. This paper is organized
as follows: Section II presents a general method for carrying out an operator
redefinition that provides an exponential reduction in the degree of coupling
of the system. In Section III this method is applied to the case of a compact
${\rm U}(1)$ lattice gauge theory with periodic boundary conditions, where a
global constraint arises from the magnetic Gauss’ law. The conclusions are
presented in Section IV.

##  II Exponential reduction in the degree of coupling

It is instructive to start with a general Hamiltonian that has a maximal
degree of coupling, in order to see how the exponential scaling arises. This
Hamiltonian is given by

$\displaystyle H$
$\displaystyle=\sum_{i=1}^{N}\left(f\left[\mathcal{Q}_{i}\right]+g\left[\mathcal{P}_{i}\right]\right)+F\left[\sum_{i=1}^{N}\mathcal{Q}_{i}\right]+G\left[\sum_{i=1}^{N}\mathcal{P}_{i}\right]$
$\displaystyle\equiv
H_{\mathcal{Q}}[\mathcal{Q}_{1},...,\mathcal{Q}_{N}]+H_{\mathcal{P}}[\mathcal{P}_{1},...,\mathcal{P}_{N}]\,,$
(2)

where $f,g,F,$ and $G$ are arbitrary functions; note that this manuscript will
always use square brackets to denote arguments of these functions, for ease of
legibility. Additionally, $\mathcal{Q}_{i}$ and $\mathcal{P}_{i}$ are
operators acting on single lattice sites, and $N$ is the system size; In
lattice field theories, $N$ would be the lattice volume.

As will be explained shortly, the naive computational cost of implementing a
given Hamiltonian scales exponentially with its DoC. The DoC for the first two
terms, which only depend on a single local operator
$\mathcal{O}_{i}\in\\{\mathcal{Q}_{1},...,\mathcal{Q}_{N},\mathcal{P}_{1},...,\mathcal{P}_{N}\\}$,
is

$\displaystyle{\rm DoC}[f;1]={\rm DoC}[g;1]=1\,,$ (3)

where we use the notation set in Eq. (1) that $\text{DoC}[X;n]$ is the Degree
of Coupling for the function $X$, which depends on $n$ operators. This is not
true for the last two terms, where the DoC depends on the number of operators
in the system $N$ and whether either $F$ or $G$ are non-polynomial. In
particular, the DoC for the third and fourth term is given by

$\displaystyle{\rm DoC}[F;N]$ $\displaystyle=\min\left[N,\deg(F)\right],$
$\displaystyle{\rm DoC}[G;N]$ $\displaystyle=\min\left[N,\deg(G)\right]\,,$
(4)

respectively. As it is the DoC of the total Hamiltonian that determines the
complexity of implementation, the relevant DoC is given by the maximum of
these two and so

$\displaystyle{\rm
DoC}=\min\left[N,\max\left[\deg(F),\deg(G)\right]\right]\,.$ (5)

Therefore, as long as either $F$ or $G$ is non-polynomial, the DoC of the
total Hamiltonian is $N$ and, as will be shown, the computational cost of
implementing this Hamiltonian will scale exponentially with the system size
$N$. However, this analysis also points towards a method for reducing this
scaling. Namely, there is a maximal imbalance between the DoC of the different
terms in the Hamiltonian as it is currently written. However, a change of
basis that better balances the DoC between the terms with ${{\rm DoC}=1}$ and
${{\rm DoC}=N}$ such that the maximum DoC is reduced would also reduce the
overall cost of implementing $H$. The main focus of this work is to present a
particular operator basis change such that the DoC of any given term in the
Hamiltonian is $\mathcal{O}(\log_{2}N)$, which reduces the computational cost
of implementing $H$ from exponential in $N$ to polynomial in $N$.

The naive computational cost of implementing this Hamiltonian can be estimated
by specifying the properties of the operators $\mathcal{Q}_{i}$ and
$\mathcal{P}_{i}$. In particular, it is assumed that $\mathcal{Q}$ and
$\mathcal{P}$ are conjugate variables, such that they obey the commutation
relations

$\displaystyle\left[\mathcal{Q}_{i},\mathcal{P}_{j}\right]=i\,\delta_{ij}\,.$
(6)

This implies that the basis $\ket{\mathcal{Q}_{i}}$ and the basis
$\ket{\mathcal{P}_{i}}$ are related by Fourier transform and that operators on
different sites commute. Therefore, one method for simulating the Trotterized
time-evolution of this Hamiltonian on a quantum computer is to implement the
matrices

$\displaystyle\innerproduct{\mathcal{Q}_{i}|e^{iH_{\mathcal{Q}}}}{\mathcal{Q}_{i}|e^{iH_{\mathcal{Q}}}}{\mathcal{Q}_{j}}\,\,\text{and}\,\,\innerproduct{\mathcal{P}_{i}|e^{iH_{\mathcal{P}}}}{\mathcal{P}_{i}|e^{iH_{\mathcal{P}}}}{\mathcal{P}_{j}}\,,$
(7)

and use a (fast) Fourier Transform to switch between the two bases.

Assuming that the each operator $\mathcal{Q}_{i}$ and $\mathcal{P}_{j}$ acts
on a Hilbert space spanned by $n_{q}$ qubits (at a given lattice site), the
dimension of the spanned Hilbert space is of size $2^{n_{q}}$. It is known
that the number of gates required to implement an arbitrary diagonal matrix of
dimension $2^{n_{q}}$ is exponential in $n_{q}$; for example, without using
ancillary qubits, the count is $2^{n_{q}+1}-3$ gates [28]. Therefore,
exponentiation of each of the terms with ${\rm DoC}=1$ in Eq. (2) requires
$\mathcal{O}(2^{n_{q}})$ gates. Implementing all of the terms with ${\rm
DoC}=1$ requires $\mathcal{O}(N2^{n_{q}})$ gates; however, as these gates can
be run in parallel due to the assumption in Eq. (6), the depth of the circuit
is therefore $\mathcal{O}(2^{n_{q}})$.

Unfortunately, this is no longer true for the terms with maximal DoC, as they
tie together all $N$ operators and thus require all $n_{q}N$ qubits.
Therefore, implementing the term with maximal DoC requires
$\mathcal{O}(2^{n_{q}N})$ gates. This is the origin of the exponential scaling
with $N$ that makes implementing this constrained Hamiltonian so costly.

###  II.1 General reduction method

The exponential scaling in system size, resulting from a term that is a
function of all the operators, can be broken by a judiciously chosen change of
operator basis. One can define a new set of operators, via

$\displaystyle\mathcal{Q}^{\prime}_{i}=\mathcal{W}_{ij}\mathcal{Q}_{j}\qquad\mathcal{P}^{\prime}_{i}=\mathcal{W}_{ij}\mathcal{P}_{j}\,,$
(8)

where $\mathcal{W}$ is an orthogonal matrix. These new set of operators obey
the same commutation relations

$\displaystyle\left[\mathcal{Q}^{\prime}_{i},\mathcal{P}^{\prime}_{j}\right]$
$\displaystyle=\mathcal{W}_{ik}\mathcal{W}_{jl}[\mathcal{Q}_{k},\mathcal{P}_{l}]$
$\displaystyle=i\mathcal{W}_{ik}\mathcal{W}_{jl}\delta_{kl}$
$\displaystyle=i\mathcal{W}_{ik}\mathcal{W}_{jl}$
$\displaystyle=i\delta_{ij}\,,$ (9)

where the last line is a consequence of the orthogonality of the
$\mathcal{W}$. This implies that the eigenbases of $\mathcal{Q}^{\prime}_{i}$
and $\mathcal{P}^{\prime}_{i}$ are still related by a Fourier transform. If
the matrix $\mathcal{W}$ is chosen such that not all of the operators appear
within the same function $F$ or $G$, then the DoC of these terms will be
reduced. However, this will also increase the DoC of the terms involving the
functions $f$ and $g$. Therefore, $\mathcal{W}$ has to be judiciously chosen
to balance these two effects. Roughly speaking, in order to break the
exponential scaling, the maximum DoC for any one term must scale
logarithmically in $N$ (or as a polynomial of $\log_{2}N$). The remaining
discussion will demonstrate how this can be achieved. In particular, in this
subsection and the following, a specific form for $\mathcal{W}$ will be
proposed that will then be proven to provide this necessary tuning.

Consider the orthogonal matrix $\mathcal{W}$ to be block diagonal in structure

$\displaystyle\mathcal{W}=\begin{pmatrix}\Omega^{(1)}&0&0&\dots&0\\\
0&\Omega^{(2)}&0&\dots&0\\\ \vdots&\dots&\ddots&\vdots&\vdots\\\
0&0&0&\dots&\Omega^{(N_{S})}\\\ \end{pmatrix}\,,$ (10)

where $N_{S}$ is the number of sub-blocks and $\Omega^{(i)}$ are orthogonal
matrices of dimension $d_{(i)}$. With this, the Hamiltonian $H_{\mathcal{Q}}$
in the rotated basis (indicated by the superscript “rot”) can be written as

$\displaystyle H^{\text{rot}}_{\mathcal{Q}}=$
$\displaystyle\sum_{i=1}^{N}f\left[\sum_{j=1}^{N}\mathcal{W}_{ij}\mathcal{Q}_{j}\right]+F\left[\sum_{i=1}^{N}\sum_{j=1}^{N}\mathcal{W}_{ij}\mathcal{Q}_{j}\right]$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{N_{S}}\sum_{k=1}^{d_{(i)}}f\left[\sum_{j=1}^{d_{(i)}}\Omega^{(i)}_{kj}\mathcal{Q}_{D_{(i)}+j-1}\right]$
$\displaystyle+F\left[\sum_{i=1}^{N_{s}}\sum_{k,j=1}^{d_{(i)}}\Omega^{(i)}_{kj}\mathcal{Q}_{D_{(i)}+j-1}\right]\,,$
(11)

where $D_{(i)}$ labels the row and/or column in the matrix $\mathcal{W}$ where
a sub-block $\Omega^{(i)}$ begins,

$\displaystyle D_{(i)}=1+\sum_{n=1}^{i-1}d_{(n)}\,,$ (12)

and the summation over repeated indices has been made explicit for clarity.
The Hamiltonian $H_{\mathcal{P}}^{\text{rot}}$ will take on a corresponding
form. In order to reduce the maximal DoC, the rotation matrix needs to reduce
the number of terms within the function $F$. This can be achieved by choosing
all elements of the first column of each $\Omega^{(i)}$ to be identical; the
veracity of this choice can be seen in the following way. Since all columns of
orthogonal matrices are vectors of unit length, if all the elements of each
$\Omega^{(i)}$ are identical, this implies
$\Omega^{(i)}_{k1}=1/\sqrt{d_{(i)}}$. Additionally, orthogonality implies the
dot product of any two different columns is zero. Because the first column
contains only identical entries, this implies the sum of all other individual
columns must be zero, i.e.

$\displaystyle\sum_{i=1}^{d_{(k)}}\Omega^{(k)}_{ij}=\sqrt{d_{(k)}}\,\delta_{1j}\,.$
(13)

With this choice, the argument of $F$ becomes

$\displaystyle\sum_{i=1}^{N_{s}}\sum_{k,j=1}^{d_{(i)}}\Omega^{(i)}_{kj}\mathcal{Q}_{D_{(i)}+j-1}$
$\displaystyle=\sum_{i=1}^{N_{s}}\sum_{j=1}^{d_{(i)}}\sqrt{d_{(i)}}\delta_{1j}\mathcal{Q}_{D_{(i)}+j-1}$
$\displaystyle=\sum_{i=1}^{N_{s}}\sqrt{d_{(i)}}\mathcal{Q}_{D_{(i)}}\,,$ (14)

and so the term that had the maximal DoC becomes

$\displaystyle F\left[\sum_{i=1}^{N}\mathcal{Q}_{i}\right]\Rightarrow
F\left[\sum_{i=1}^{N_{S}}\sqrt{d_{(i)}}\,\mathcal{Q}_{D_{(i)}}\right]\,.$ (15)

This implies that the DoC is no longer $N$ but $N_{S}$ and therefore, the
number of gates needed to implement this term is reduced to
$\mathcal{O}(2^{n_{q}N_{S}})$. However, the exponential volume scaling will
persist unless $N_{S}\lesssim(\log_{2}N)^{p}$, for some power $p$.

While the change of basis has decreased the DoC for the function $F$, it has
done so at the cost of increasing the DoC for the function $f$. In particular,
the $i^{\text{th}}$ term with $\rm{DoC}=1$ becomes

$\displaystyle
f\left[\mathcal{Q}_{i}\right]\Rightarrow\sum_{k=1}^{d_{(i)}}f\left[\sum_{j=1}^{d_{(i)}}\Omega^{(i)}_{kj}\mathcal{Q}_{D_{(i)}+j-1}\right]\,,$
(16)

and therefore, the DoC of the $i^{\text{th}}$ term is now given by the number
of non-zero entries in the $j^{\text{th}}$ row of $\Omega^{(i)}$, defined as
$\varphi^{(i)}_{j}$. This implies that the number of gates required to
implement Eq. (16) for all $N$ sites scales as

$\displaystyle\mathcal{O}\left(\sum_{i=1}^{N_{s}}\sum_{j=1}^{d_{(i)}}2^{n_{q}\varphi^{(i)}_{j}}\right)\,,$
(17)

which is exponential in $\varphi_{j}^{(i)}$.

To recap, in the original operator basis, the DoC of $H_{\mathcal{Q}}$ is
given by

$\displaystyle\rm DoC$ $\displaystyle=\max\left[\rm DoC\left[f;1\right],\rm
DoC\left[F;N\right]\right]$ $\displaystyle=\min\left[N,\deg(F)\right],$ (18)

while in the rotated basis it is

DoC
$\displaystyle=\max\left[\text{DoC}\left[f;\max(\varphi_{j}^{(i)})\right],\text{DoC}\left[F;N_{S}\right]\right],$
(19)

where $\varphi_{j}^{(i)}$ is a measure of the sparcity of $\Omega^{(i)}$.
Considering the case where both $f$ and $F$ are non-polynomial functions and
thus have $\deg(f)=\deg(F)=\infty$, the DoC in the rotated basis is

DoC $\displaystyle=\max\left[\max(\varphi_{j}^{(i)}),N_{S}\right].$ (20)

This implies that the exponential volume scaling can be reduced to polynomial
if one finds a basis such that
$\max\left[\max(\varphi_{j}^{(i)}),N_{S}\right]=\mathcal{O}(\log_{2}N)^{p}$
for some power $p$; the same is true for $H_{\mathcal{P}}$ as it has the same
structure as $H_{\mathcal{Q}}$. The following section will prove that it is
always possible to achieve this tuning.

###  II.2 Weaved rotation matrices

This section will discuss how to construct rotation matrices which satisfy the
criteria

$\displaystyle\varphi_{j}\leq\lceil\log_{2}d\rceil+1\,.$ (21)

Throughout this section it is generally assumed that each sub-block of
$\mathcal{W}$ is of similar size, $d_{(i)}\sim d$, and also that every row of
each sub-block has a comparable number of non-zero elements,
$\varphi^{(i)}_{j}\sim\varphi$. These assumptions, combined with the above
criteria, imply that the number of gates needed to implement Eq. (11) is

$\displaystyle\mathcal{O}\left(N^{n_{q}}+N\left(\frac{N}{\log_{2}N}\right)^{n_{q}}\right)\,,$
(22)

which scales polynomially, not exponentially, with the system size. As
$H_{\mathcal{P}}$ has the same structure, it will have the same reduction in
the number of gates needed to implement its complex exponential. Note that the
exact choice of $N_{S}$, i.e. whether $\log_{2}N$ is rounded up or down, will
have $\mathcal{O}(1)$ effects on the number of gates; this can be readily
optimized.

Let the matrix $\Omega^{(i)}_{M}$ be any orthogonal matrix of dimension $M$
whose first column has all entries equal to $1/\sqrt{M}$; the index $(i)$ is
used to denote a specific matrix as $\Omega_{M}$ of size $M$ is not unique.
Some examples of such matrices are

$\displaystyle\Omega^{(a)}_{4}$
$\displaystyle=\begin{pmatrix}\frac{1}{2}&-\frac{1}{\sqrt{2}}&-\frac{1}{2}&0\\\
\frac{1}{2}&\frac{1}{\sqrt{2}}&-\frac{1}{2}&0\\\
\frac{1}{2}&0&\frac{1}{2}&-\frac{1}{\sqrt{2}}\\\
\frac{1}{2}&0&\frac{1}{2}&\frac{1}{\sqrt{2}}\end{pmatrix}$
$\displaystyle\Omega^{(b)}_{4}$
$\displaystyle=\frac{1}{2}\begin{pmatrix}1&1&1&1\\\ 1&1&-1&-1\\\ 1&-1&-1&1\\\
1&-1&1&-1\end{pmatrix}\,.$ (23)

This section proves that there is a subset of the matrices $\Omega_{M}$ that
obey

$\displaystyle\eta(\Omega_{M})=\lceil\log_{2}(M)\rceil+1\quad\forall\quad
M\in\mathbb{Z}^{+}\,.$ (24)

Here $\eta(\mathbb{M})$ is the sparsity of the matrix $\mathbb{M}$, defined as
the largest number of non-zero entries that appear in any individual row of
the matrix $\mathbb{M}$, i.e.
$\eta(\mathbb{M})=\text{max}_{j}(\varphi_{j}(\mathbb{M}))$, where
$\varphi_{j}(\mathbb{M})$ was previously defined as the number of non-zero
entries in the $j^{\text{th}}$ row of a matrix $\mathbb{M}$. For example, in
the case above,

$\displaystyle\eta(\Omega_{4}^{(a)})=3\qquad\eta(\Omega_{4}^{(b)})=4\,.$ (25)

With this definition, $\Omega^{(a)}_{4}$ obeys Eq. (24) while
$\Omega^{(b)}_{4}$ does not. Matrices that obey Eq. (24) will be called weaved
matrices, $W_{M}$, where the subscript again labels the dimension of the
square orthogonal matrix. The reason behind this naming convention will be
made clear when their explicit construction is discussed. The proof will
proceed in two parts. The first focuses on constructing weaved matrices of
dimensions $M=2^{m},m\in\mathbb{Z}^{0+}$, while the second focuses on
constructing matrices of all other dimensions.

For the first part, let $T^{(i,j)}_{M}(\theta)$ be a rotation matrix that
performs a 2-dimensional rotation by an angle $\theta$ in an $M$-dimensional
vector space in the plane spanned by the $i^{\text{th}}$ and $j^{\text{th}}$
component of the vector space. In other words, one can write

$\displaystyle T^{(i,j)}_{M}(\theta)$
$\displaystyle=\begin{pmatrix}\mathbbm{1}_{i-1}&0&0&0&0\\\
0&\cos\theta&0&-\sin\theta&0\\\ 0&0&\mathbbm{1}_{j-i-1}&0&0\\\
0&\sin\theta&0&\cos\theta&0\\\ 0&0&0&0&\mathbbm{1}_{M-j}\\\ \end{pmatrix}\,,$
(26)

where $\mathbbm{1}_{i}$ is the $i\times i$ identity matrix. A few relevant
examples are

$\displaystyle T^{(1,2)}_{2}(\theta)$
$\displaystyle=\begin{pmatrix}\cos\theta&-\sin\theta\\\
\sin\theta&\cos\theta\end{pmatrix}\,,$ $\displaystyle T^{(1,3)}_{4}(\theta)$
$\displaystyle=\begin{pmatrix}\cos\theta&0&-\sin\theta&0\\\ 0&1&0&0\\\
\sin\theta&0&\cos\theta&0\\\ 0&0&0&1\end{pmatrix}\,.$ (27)

This matrix is used to create $W_{2^{m}}$, assuming that $W_{2^{m-1}}$ has
already been constructed. In particular, $T$ weaves together two copies of
$W_{2^{m-1}}$ via

$\displaystyle W_{2^{m}}=\begin{pmatrix}W_{2^{m-1}}&0\\\
0&W_{2^{m-1}}\end{pmatrix}\times
T^{(1,1+2^{m-1})}_{2^{m}}\left(\pi/4\right)\,,$ (28)

Notice that due to the form of the matrix $T$, explicitly shown in Eq. (26),
its effect is to introduce a single non-zero entry into each row of the off-
diagonal sub-blocks. For example,

$\displaystyle
W_{4}=\begin{pmatrix}\frac{W_{2}^{(11)}}{\sqrt{2}}&W_{2}^{(12)}&-\frac{W_{2}^{(11)}}{\sqrt{2}}&0\\\
\frac{W_{2}^{(21)}}{\sqrt{2}}&W_{2}^{(22)}&-\frac{W_{2}^{(21)}}{\sqrt{2}}&0\\\
\frac{W_{2}^{(11)}}{\sqrt{2}}&0&\frac{W_{2}^{(11)}}{\sqrt{2}}&W_{2}^{(12)}\\\
\frac{W_{2}^{(21)}}{\sqrt{2}}&0&\frac{W_{2}^{(21)}}{\sqrt{2}}&W_{2}^{(22)}\end{pmatrix}$
(29)

Therefore, the $\eta$ parameter of the matrix $W_{2^{m}}$ is given by

$\displaystyle\eta(W_{2^{m}})=\eta(W_{2^{m-1}})+1\,.$ (30)

However, $W_{2^{m-1}}$ can again be created by weaving together two copies of
$W_{2^{m-2}}$ using an analogous $T$ matrix, with $m\rightarrow m-1$. By
iteration one therefore finds

$\displaystyle\eta(W_{2^{m}})$ $\displaystyle=\eta(W_{2^{0}})+m\,.$ (31)

Since $W_{2^{0}}=1$, it follows that

$\displaystyle\eta(W_{2^{m}})=1+m\,.$ (32)

Therefore $W_{2^{m}}$ satisfies Eq. (24) for any non-zero integer $m$. The
orthogonality of $W_{2^{m}}$ can also be proven inductively, as $W_{2^{m}}$ is
orthogonal as long as $W_{2^{m-1}}$ is due to the matrix $T$ being a rotation
matrix. The first few weaved matrices are given explicitly by

$\displaystyle W_{1}$ $\displaystyle=1,$ $\displaystyle W_{2}$
$\displaystyle=\frac{1}{\sqrt{2}}\begin{pmatrix}1&-1\\\ 1&1\end{pmatrix},$
$\displaystyle W_{4}$
$\displaystyle=\begin{pmatrix}\frac{1}{2}&-\frac{1}{\sqrt{2}}&-\frac{1}{2}&0\\\
\frac{1}{2}&\frac{1}{\sqrt{2}}&-\frac{1}{2}&0\\\
\frac{1}{2}&0&\frac{1}{2}&-\frac{1}{\sqrt{2}}\\\
\frac{1}{2}&0&\frac{1}{2}&\frac{1}{\sqrt{2}}\end{pmatrix}.$ (33)

To summarize, the results presented so far show that it is possible to
construct an orthogonal matrix of dimension $2^{m}$, with $m$ any positive
integer, whose first column has all entries equal to $1/\sqrt{2^{m}}$ and

$\displaystyle\eta(W_{2^{m}})=m+1\,.$ (34)

The second part of the proof shows how to construct weaved matrices for
general values of $M\neq 2^{m}$. This makes use of the fact that any integer
can be written in a binary form with $K=\lceil\log_{2}(M)\rceil$ bits, such
that one can define a set $\tilde{m}(M)$ containing the position of the non-
zero bits. For example, $\tilde{m}(5)=\\{0,2\\}$, while
$\tilde{m}(10)=\\{1,3\\}$. Furthermore, $\tilde{m}_{i}$ denotes the
$i^{\rm{th}}$ entry in $\tilde{m}(M)$. With this, $W_{M}$, for $M\neq 2^{m}$,
is defined to be

$\displaystyle W_{M}$
$\displaystyle=\begin{pmatrix}W_{2^{\tilde{m}_{1}}}&0&0\\\ 0&\ddots&0\\\
0&0&W_{2^{\tilde{m}_{k}}}\end{pmatrix}\times\prod_{j=1}^{k-1}T^{\left(1,1+b_{j}(\tilde{m})\right)}_{M}(\Theta_{j})\,,$
(35)

where $k\equiv\dim\tilde{m}(M)$ and

$\displaystyle b_{j}(\tilde{m})=\sum_{i=1}^{j}2^{\tilde{m}_{i}}\,,$ (36)

which is the binary representation that only keeps the first $j$ digits. The
matrix multiplication of the product reads left to right, $T^{(i,j)}_{M}$ is
the same rotation matrix defined above and $\theta_{j}$ is given by

$\displaystyle\theta_{j}=\arccos\left[\left(1+\left[2^{\tilde{m}_{j+1}}/b_{j}(\tilde{m})\right]\right)^{-1/2}\right]\,.$
(37)

Note that the product of $T$ matrices can be simplified to

$\displaystyle\prod_{j=1}^{k-1}T^{\left(1,1+b_{j}(\tilde{m})\right)}_{M}(\Theta_{j})=$
(38)
$\displaystyle\begin{pmatrix}\scriptstyle\prod_{j=1}^{k-1}\cos(\Theta_{j})&0&-\sin\Theta_{1}&0&\cdots&-\sin\Theta_{k-1}&0\\\
0&\mathbbm{1}_{\Delta_{1}}&0&0&\cdots&0&0\\\
\sin\Theta_{1}&0&\cos\Theta_{1}&0&\cdots&0&0\\\
0&0&0&\mathbbm{1}_{\Delta_{2}}&\cdots&0&0\\\
\vdots&\vdots&\vdots&\vdots&\ddots&0&0\\\
\sin\Theta_{k-1}&0&0&0&\cdots&\cos\Theta_{k}&0\\\
0&0&0&0&\cdots&0&\mathbbm{1}_{2^{\tilde{m}_{k}}-1}\end{pmatrix}\,,$

where $\Delta_{i}$ measures the difference between the dimension of
‘neighboring’ $W_{2^{\tilde{m}_{i}}}$,

$\displaystyle\Delta_{i}\equiv 2^{\tilde{m}_{i+1}}-2^{\tilde{m}_{i}}-1\,.$
(39)

Note that because $W_{M}$ is constructed from products of orthogonal matrices,
$W_{M}$ is itself orthogonal.

It will now be shown that $W_{M}$ obeys the sparsity constraint

$\displaystyle\eta(W_{M})=\lceil\log_{2}(M)\rceil+1\,.$ (40)

As shown in the first part of the proof, the sparsity of each of the
$W_{2^{\tilde{m}_{i}}}$ sub-blocks is given by $1+\tilde{m}_{i}$. Each of
these blocks is multiplied non-trivially by at most
$\lfloor\log_{2}M\rfloor-\tilde{m}_{i}+1$ $T$-matrices from the right, each of
which increases the sparsity by one. Furthermore, the largest $\tilde{m}_{i}$
is given by

$\displaystyle\tilde{m}_{k}=\lfloor\log_{2}M\rfloor\,,$ (41)

and its corresponding $W_{2^{\tilde{m}_{k}}}$ is multiplied non-trivially by
exactly one $T$-matrix. Therefore

$\displaystyle\eta(W_{M})$ $\displaystyle=1+\left(\tilde{m}_{k}+1\right)$
$\displaystyle=1+\lceil\log_{2}M\rceil\,,$ (42)

for any value of $M$.

To summarize, it is always possible to construct an orthogonal matrix of
dimension $M$ whose first column has all entries equal to $1/\sqrt{M}$ and
whose rows contain at most $\lceil\log_{2}M\rceil+1$ non-zero entries. It is
instructive to present a visualization of the dramatic change in the DoC if
the procedure described in the previous section is carried out using weaved
matrices. In Fig. 1, the DoC for each operator in a system that has a total of
sixteen operators is presented. Specifically, Fig. 1a shows the coupling
between operators in the original basis. Note that this is a maximally coupled
system. On the other hand, Fig. 1b shows the coupling between sites in the
weaved basis, which has a significant reduction.

(a)  Degree of coupling in the original basis. The system has maximum DoC,
resulting in circuit depths that are exponential in the volume.

(b) Degree of connectivity in the weaved basis. This reduction in connectivity
results in a circuit depth that scales polynomially in the volume.

Figure 1: Diagrammatic depiction of the Degree of Connectivity for a system
with 16 total operators. Fig. 1a shows the result for a Hamiltonian with the
original constraint giving rise to maximal ${\rm DoC}=16$. Fig. 1b illustrates
the DoC after the operator basis has been transformed to reduce the degree of
connectivity, using the weaved operation. The degree of connectivity is now
$\rm DoC=4$, coming from the terms involving $\mathcal{O}_{1}$,
$\mathcal{O}_{5}$, $\mathcal{O}_{9}$ and $\mathcal{O}_{13}$. Note that both
Hamiltonians have the same spectrum, up to digitization effects. The operator
$\mathcal{O}_{i}$ could be $\mathcal{P}_{i}$ or $\mathcal{Q}_{i}$.

###  II.3 Classical Computational Cost of Constructing Weaved Matrices

In this section, a derivation of the classical computation cost required to
construct weaved matrices is given. It will show that the classical cost
required to construct the weaved matrix $W_{M}$ is $\mathcal{O}(M\log_{2}M)$.

Recall that the naive cost of multiplying together two arbitrary square
matrices of dimension $M$ is $\mathcal{O}(M^{3})$. This can be seen by noting
that each of the $\mathcal{O}(M^{2})$ entries in the new matrix requires
$\mathcal{O}(M)$ calculations. Note that there are advanced algorithms that
reduce the scaling to less than cubic [29, 30]; however, this section will
assume that the computational cost scales like the total number of simple
algebraic manipulations.

The computational cost of multiplying matrices is reduced if the matrices are
sparse, as the zero entries can simply be ignored. Therefore, if the matrices
are exceptionally sparse, the cost can be dramatically reduced. The cost of
constructing a weaved matrix of dimension $W_{2^{m}}$ with $m$ a positive
integer greater than 1 is derived first. By construction this matrix has only
$2^{m}(m+1)$ non-zero entries and furthermore, the first column has all
$2^{m}$ entries non-zero and therefore the remainder of the matrix contains
only $2^{m}m$ non-zero entries. To estimate the computational cost of
constructing this matrix, recall that the construction of $W_{2^{m}}$ is done
by weaving together two copies of $W_{2^{m-1}}$, as given in Eq. (28). This
means that the $2^{m-1}(m-1)$ non-zero entries in each copy of $W_{2^{m-1}}$,
excluding the first column, are multiplied by the identity matrix, requiring
$2^{m}(m-1)$ algebraic multiplications. Next, each first column of
$W_{2^{m-1}}$ is multiplied by $\cos\pi/4$ and $\pm\sin\pi/4$, where the sign
on the $\sin\pi/4$ depends on whether it is the lower or upper $W_{2^{m-1}}$.
This requires $2\times 2\times 2^{m-1}$ algebraic multiplications. Lastly, due
to the structure of $T$ shown in Eq. (26), these two procedures factorize and
so there is no addition that has to happen (unlike when multiplying generic
matrices). Therefore, the total number of algebraic manipulations to construct
$W_{2^{m}}$, assuming that $W_{2^{m-1}}$ is already constructed, is
$2^{m}(m+1)$, which is equal to the total number of non-zero entries in
$W_{2^{m}}$. If $W_{2^{m-1}}$ is not constructed, then it will take $2^{m-1}m$
algebraic manipulations. To construct $W_{2^{m}}$, starting from $W_{1}$ thus
requires

$\displaystyle\sum_{j=1}^{m}2^{j}(j+1)=m2^{m+1}$ (43)

algebraic manipulations. The total classical computational cost of
constructing $W_{2^{m}}$ is therefore $\mathcal{O}(m2^{m})$.

Similar arguments can be used to understand the scaling of constructing
$W_{M}$, for $M$ an arbitrary positive integer. Constructing all of the
necessary $W_{2^{\tilde{m}_{i}}}$ will take $\mathcal{O}(M\log_{2}M)$ as the
maximum $\tilde{m}$ is given by $\lfloor\log_{2}M\rfloor$. The construction of
$W_{M}$ is again constructed by multiplying with the product of $T$ matrices
shown in Eq. (II.2). Repeating similar steps as before still results in the
total computational cost of computing $W_{M}$ scaling as
$\mathcal{O}(M\log_{2}M)$.

###  II.4 Comments on generalizations of this approach

Note that the discussion so far assumed that the two operators obeyed the
canonical commutation relations. This fact was important to ensure that the
bases of $\mathcal{Q}^{\prime}_{i}$ and $\mathcal{P}^{\prime}_{i}$ are still
related by a Fourier transform, and was not material in reducing the DoC
between the operators. One therefore expects that similar ‘change of basis’
methods would also work for Hamiltonians which are constructed out of
operators that do not satisfy canonical commutation relations; this would
require finding an efficient implementation of the new relation between the
two basis $\mathcal{Q}^{\prime}_{i}$ and $\mathcal{P}^{\prime}_{i}$.

Besides this, there are a few important caveats about the applicability and
usefulness of this method applied to these other cases. The first, most
important caveat is that, as currently written, this method is only exact when
applied to operators that have an infinite dimensional spectrum. This can be
understood by realizing that when performing Kronecker sums of operators, the
dimensionality of the space that they span is increased. However, if this
space is already infinite, it does not make a difference. The generalization
and application to discrete operators has interesting applications and will be
carried out in future work.

Second, for polynomial functions $F$ and $G$ of order $p$, no term contains
more than $p$ operators multiplied together. Thus, the DoC is in fact no
larger than the order of the polynomial $p$, and the number of gates
required222The case with $p=1$ and $p=2$ was shown in Ref. [31]; the
generalization to arbitrary $p$ can be understood by counting the number of
individual terms in $\bigotimes_{i=1}^{p}x_{i}$, with
$x=\frac{x_{\text{max}}}{2^{n_{q}}-1}\sum_{j=1}^{n_{q}}2^{j}\sigma^{z}_{j}$.
is $\mathcal{O}(n_{q}^{p})$. However, if $p$ is large, it still might be
beneficial to utilize the change of basis in order to decrease the coefficient
in front of the polynomial scaling.

##  III Application to Compact U(1) Gauge Theory

An example of where this change of basis is of fundamental importance is for
an implementation of the full gauge-fixed compact 2+1 dimensional pure U(1)
gauge theory. In this section, after briefly reviewing the formulation, the
change of basis is performed and the dramatic reduction in the number of gates
needed to implement time evolution of this Hamiltonian using Suzuki-Trotter
methods is discussed.

###  III.1 Review of lattice formulation of compact U(1) gauge theory

The Hamiltonian of a ${\rm U}(1)$ gauge theory in two spatial dimensions is
given by

$\displaystyle H=\int{\rm d}^{2}x\left(\vec{E}(x)^{2}+B(x)^{2}\right)\,,$ (44)

where the electric and magnetic fields $\vec{E}(x)$ and $B(x)$, respectively,
are related to the vector potential via

$\displaystyle\vec{E}(x)$ $\displaystyle=\partial_{t}\vec{A}(x)$
$\displaystyle B(x)$ $\displaystyle=\vec{\nabla}\times\vec{A}(x)\,.$ (45)

Note that in two dimensions the curl of a vector field yields a pseudoscalar,
while in the more familiar three dimensions it yields a vector. In a theory
without matter, as considered here, Gauss’ law gives a constraint on the
electric field $\vec{\nabla}\cdot\vec{E}=0$.

Putting the theory onto a spatial lattice, the Hamiltonian can still be split
into a contribution from electric and magnetic fields

$\displaystyle H=H_{E}+H_{B}\,,$ (46)

where electric fields are defined on the links that connect the lattice sites,
and magnetic fields are defined on plaquettes, which are constructed out of
four links. The Hamiltonian can be written in a gauge-redundancy free manner
by going into the dual basis formulation [19, 20, 21, 22]. In particular, the
electric links are exchanged for rotor variables, which are plaquette
variables, via

$\displaystyle\vec{E}_{\ell}=\vec{\nabla}_{L}\times R_{p}\,,$ (47)

where $\vec{\nabla}_{L}$ is the differential operator, defined on a lattice.
The magnetic field is left unaltered, defined on the plaquettes as $B_{p}$. In
a quantum theory, $R_{p}$ and $B_{p}$ are promoted to operators, which satisfy
the canonical commutation relations

$\displaystyle\left[\hat{B}_{p},\hat{R}_{p^{\prime}}\right]=i\,\delta_{pp^{\prime}}\,.$
(48)

In a pure gauge theory with no charges, the (local) Gauss law constraint

$\displaystyle\vec{\nabla}_{L}\cdot\vec{\hat{E}}_{\ell}=0$ (49)

is automatically satisfied. However, for periodic boundary conditions, there
is a global constraint that can be thought of as a magnetic Gauss law, stating
that the total magnetic field through a closed surface vanishes

$\displaystyle\sum_{p}\hat{B}_{p}=0\,.$ (50)

Putting this information together, the electric Hamiltonian is given by

$\displaystyle\hat{H}_{E}=\frac{g^{2}}{2a}\sum_{p=1}^{N_{p}}(\vec{\nabla}_{L}\times\hat{R}_{p})\,.$
(51)

where $g$ is the dimensionless gauge coupling; note that this Hamiltonian is
bilinear and also nearest neighbor in the rotor fields. Note that $N_{p}$ is
the number of independent plaquettes, given by $N_{p}=N_{x}N_{y}-1$ where
$N_{x(y)}$ is the number of lattice sites in the $\hat{x}(\hat{y})$ direction.
The magnetic Hamiltonian, for the compact version of the theory is given by

$\displaystyle\hat{H}_{B}=-\frac{1}{2a\,g^{2}}\left(\sum_{p=1}^{N_{p}}\cos\left[\hat{B}_{p}\right]+\cos\left[\sum_{p=1}^{N_{p}}\hat{B}_{p}\right]\right)\,,$
(52)

where an overall constant shift has been dropped.

The implementation of this Hamiltonian onto a digital quantum computer
requires the introduction of a digitization and truncation scheme. This was
done in Ref. [32], but the exact details are not necessary for carrying out
the change of basis procedure described in Sec. II.1. The only relevant detail
is that in this scheme, the operators $\hat{B}_{p}$ and $\hat{R}_{p}$ have
equally spaced eigenvalues and their eigenbases are related by Fourier
transform. Additionally, since the functions $f$ and $F$ are periodic, care
needs to be taken that eigenvalues of the new operators do not exceed the
range $-\pi$ to $\pi$.

Before implementing the change of basis procedure, it is again important to
highlight that there are many other proposed Hamiltonian formulations and
computational strategies for simulating U(1) (or other gauge theories) on a
quantum computer [20, 33, 21, 22, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44,
45, 46, 47, 48, 49, 50, 51, 52]. Not all of these proposals depend on a fully
gauge-fixed Hamiltonian. Instead, they utilize other methods for ensuring that
Gauss’ Laws are ultimately satisfied. It is currently an open question of what
approach is best, and in fact, the answer may change as quantum computers
evolve beyond the Noisy-Intermediate Scale Quantum (NISQ)-era [53]. The
purpose of this section is to demonstrate how one of the main drawbacks of
using fully gauge-fixed Hamiltonians, the high DoC, may be ameliorated.

###  III.2 Utilizing the Change of Basis Procedure

The Hamiltonian for the compact ${\rm U}(1)$ theory presented above is of the
general form given in Eq. (2) with

$\displaystyle f(x)=F(x)=\cos(x)\qquad g(x)=G(x)=x^{2}\,.$ (53)

As $g$ and $G$ are both quadratic in the operators, the number of gates needed
to implement their time evolution scales quadratically with $n_{q}$, the
number of qubits used to represent each operator. Because there are $N_{p}$
such terms, the number of gates required to implement $g$ and $G$ is
$\mathcal{O}(N_{p}n_{q}^{2})$. However, as both $f$ and $F$ are non-polynomial
functions, the gates needed to implement their time evolution will naively
scale exponentially with the number of qubits each operator acts on. For all
$N_{p}$ of the $f$ terms the scaling would then be $\mathcal{O}(2^{n_{q}})$,
leading to a total scaling of $\mathcal{O}(N_{p}2^{n_{q}})$. The function $F$
acts on the entire system of qubits, and therefore scales as
$\mathcal{O}(2^{n_{q}N_{p}})$, which is exponential in $N_{p}$. However, this
naive scaling can be reduced using the methods developed in Sec. II, leading
to a dramatic reduction of the number of gates needed.

Before providing some numerical examples of the improvement in scaling, note
that carrying out a change of basis in the magnetic Hamiltonian will alter the
electric Hamiltonian. Because $H_{E}$ consists of only bilinear terms,
performing a basis change can introduce up to a maximum of $N_{p}^{2}$ such
terms, leading to a worst case scaling of $\mathcal{O}(N_{p}^{2}n_{q}^{2})$ in
the weaved basis. However, a potential growth in cost in the electric
Hamiltonian will be more than offset by the breaking of the exponential
scaling in the magnetic Hamiltonian.

To get a sense of the full scope of this improvement, one can estimate the
number of gates necessary for implementing the original versus weaved
Hamiltonians. Using $\mathcal{O}(2^{n_{q}N})$ scaling, the original operator
basis requires $\mathcal{O}(10^{5})$ gates for a $3\times 3$ lattice and
$\mathcal{O}(10^{9})$ gates for a $4\times 4$ lattice, assuming two qubits per
operator. In the weaved basis, using Eq. (22), both of these lattices require
only $\mathcal{O}(10^{2})$ gates.

##  IV Conclusion

This paper presented a general method to reduce the degree of coupling of a
Hamiltonian of the form

$\displaystyle H$
$\displaystyle=\sum_{i=1}^{N}\left(f\left[\mathcal{Q}_{i}\right]+g\left[\mathcal{P}_{i}\right]\right)+F\left[\sum_{i=1}^{N}\mathcal{Q}_{i}\right]+G\left[\sum_{i=1}^{N}\mathcal{P}_{i}\right]\,,$
(54)

which arises in systems where the Hamiltonian is given by a sum over functions
of individual operators, with a global constraint on the sum of these
operators. A well known example of such a system is a compact ${\rm U}(1)$
gauge theory that is formulated on a lattice with periodic boundary
conditions.

The method requires performing a change of operator basis, which reduces to
maximal number of terms involved in the arguments of a non-polynomial function
from $N$ to $\log_{2}N$. The matrices defining the operator redefinition are
constructed using iterative methods by weaving together lower dimensional
matrices in a particular way. These matrices can be constructed very
efficiently, with the computational cost scaling as $N\log_{2}N$.

The presented technique breaks the computational cost required to simulate
such a theory from scaling exponentially with the system size $N$ to scaling
polynomial with the system size. The expectation is that this method will have
applications beyond the example of a ${\rm U}(1)$ theory which was the main
motivation for this work. A detailed study of the implementation of the ${\rm
U}(1)$ theory is left for future work.

###### Acknowledgements.

The authors would like to thank Hank Lamm, Natalie Klco, Jesse Stryker, Tobias
Osborne and Zhiyao Li for useful discussions and comments on the manuscript.
CWB and BN are supported by the U.S. Department of Energy (DOE), Office of
Science under contract DE-AC02-05CH11231. In particular, support comes from
Quantum Information Science Enabled Discovery (QuantISED) for High Energy
Physics (KA2401032). DMG acknowledges the support of the CERN Quantum
Technology Initiative as well as thanking the Mainz Institute for Theoretical
Physics (MITP) of the Cluster of Excellence PRISMA+ (Project ID 39083149) for
its hospitality and support during the writing of this manuscript. CK is
supported by the DOE Computational Science Graduate Fellowship under award
number DE-SC0020347. This research benefited from the resources of the Oak
Ridge Leadership Computing Facility, which is a DOE Office of Science User
Facility supported under Contract DE-AC05-00OR22725

## References

  * Dirac [1964] P. A. M. Dirac,  _Lectures on quantum mechanics_ (1964).
  * Henneaux and Teitelboim [1992] M. Henneaux and C. Teitelboim,  _Quantization of Gauge Systems_ (Princeton University Press, 1992).
  * Bojowald [2010] M. Bojowald,  _Canonical Gravity and Applications: Cosmology, Black Holes, and Quantum Gravity_ (Cambridge University Press, 2010).
  * Suzuki [1976a] M. Suzuki, Commun. Math. Phys. 51, 183 (1976a).
  * Suzuki [1976b] M. Suzuki, Progress of Theoretical Physics 56, 1454 (1976b), https://academic.oup.com/ptp/article-pdf/56/5/1454/5264429/56-5-1454.pdf .
  * Trotter [1959] H. F. Trotter, Proceedings of the American Mathematical Society 10, 545 (1959).
  * Zohar and Reznik [2011] E. Zohar and B. Reznik, Phys. Rev. Lett. 107, 275301 (2011).
  * Zohar  _et al._ [2012] E. Zohar, J. I. Cirac, and B. Reznik, Phys. Rev. Lett. 109, 125302 (2012).
  * Hauke  _et al._ [2013a] P. Hauke, D. Marcos, M. Dalmonte, and P. Zoller, Physical Review X 3 (2013a), 10.1103/physrevx.3.041018.
  * Halimeh and Hauke [2020] J. C. Halimeh and P. Hauke, Phys. Rev. Lett. 125, 030503 (2020), arXiv:2001.00024 [cond-mat.quant-gas] .
  * Lamm  _et al._ [2019] H. Lamm, S. Lawrence, and Y. Yamauchi (NuQS), Phys. Rev. D 100, 034518 (2019), arXiv:1903.08807 [hep-lat] .
  * Lamm  _et al._ [2020] H. Lamm, S. Lawrence, and Y. Yamauchi (NuQS), (2020), arXiv:2005.12688 [quant-ph] .
  * Tran  _et al._ [2021] M. C. Tran, Y. Su, D. Carney, and J. M. Taylor, P. R. X. Quantum. 2, 010323 (2021), arXiv:2006.16248 [quant-ph] .
  * Banerjee  _et al._ [2012] D. Banerjee, M. Dalmonte, M. Muller, E. Rico, P. Stebler, U. J. Wiese, and P. Zoller, Phys. Rev. Lett. 109, 175302 (2012), arXiv:1205.6366 [cond-mat.quant-gas] .
  * Kasper  _et al._ [2020] V. Kasper, T. V. Zache, F. Jendrzejewski, M. Lewenstein, and E. Zohar, (2020), arXiv:2012.08620 [quant-ph] .
  * Hauke  _et al._ [2013b] P. Hauke, D. Marcos, M. Dalmonte, and P. Zoller, Phys. Rev. X 3, 041018 (2013b).
  * Kühn  _et al._ [2014] S. Kühn, J. I. Cirac, and M.-C. Bañuls, Physical Review A 90 (2014), 10.1103/physreva.90.042305.
  * Stannigel  _et al._ [2014] K. Stannigel, P. Hauke, D. Marcos, M. Hafezi, S. Diehl, M. Dalmonte, and P. Zoller, Phys. Rev. Lett. 112, 120406 (2014), arXiv:1308.0528 [quant-ph] .
  * Drell  _et al._ [1979] S. D. Drell, H. R. Quinn, B. Svetitsky, and M. Weinstein, Phys. Rev. D 19, 619 (1979).
  * Kaplan and Stryker [2020] D. B. Kaplan and J. R. Stryker, Phys. Rev. D 102, 094515 (2020), arXiv:1806.08797 [hep-lat] .
  * Haase  _et al._ [2021] J. F. Haase, L. Dellantonio, A. Celi, D. Paulson, A. Kan, K. Jansen, and C. A. Muschik, Quantum 5, 393 (2021), arXiv:2006.14160 [quant-ph] .
  * Bender and Zohar [2020] J. Bender and E. Zohar, Phys. Rev. D 102, 114517 (2020), arXiv:2008.01349 [quant-ph] .
  * Wilson [1974] K. G. Wilson, Phys. Rev. D 10, 2445 (1974).
  * Kogut and Susskind [1975] J. Kogut and L. Susskind, Phys. Rev. D 11, 395 (1975).
  * Banks  _et al._ [1976] T. Banks, L. Susskind, and J. Kogut, Phys. Rev. D 13, 1043 (1976).
  * Polyakov [1975] A. M. Polyakov, Phys. Lett. B 59, 82 (1975).
  * Polyakov [1977] A. M. Polyakov, Nucl. Phys. B 120, 429 (1977).
  * Bullock and Markov [2003] S. S. Bullock and I. L. Markov,  (2003), 10.48550/ARXIV.QUANT-PH/0303039.
  * Strassen [1969] V. Strassen, Numerische Mathematik 13, 354 (1969).
  * Alman and Williams [2020] J. Alman and V. V. Williams, “A refined laser method and faster matrix multiplication,”  (2020).
  * Klco and Savage [2019] N. Klco and M. J. Savage, Phys. Rev. A 99, 052335 (2019).
  * Bauer and Grabowska [2021] C. W. Bauer and D. M. Grabowska,  (2021), 10.48550/ARXIV.2111.08015.
  * Unmuth-Yockey [2019] J. F. Unmuth-Yockey, Phys. Rev. D 99, 074502 (2019), arXiv:1811.05884 [hep-lat] .
  * Anishetty  _et al._ [2010] R. Anishetty, M. Mathur, and I. Raychowdhury, J. Phys. A 43, 035403 (2010), arXiv:0909.2394 [hep-lat] .
  * Raychowdhury and Stryker [2020a] I. Raychowdhury and J. R. Stryker, Phys. Rev. Res. 2, 033039 (2020a), arXiv:1812.07554 [hep-lat] .
  * Raychowdhury and Stryker [2020b] I. Raychowdhury and J. R. Stryker, Phys. Rev. D 101, 114502 (2020b), arXiv:1912.06133 [hep-lat] .
  * Anishetty  _et al._ [2009] R. Anishetty, M. Mathur, and I. Raychowdhury, J. Math. Phys. 50, 053503 (2009), arXiv:0901.0644 [math-ph] .
  * Mathur  _et al._ [2010] M. Mathur, I. Raychowdhury, and R. Anishetty, Journal of Mathematical Physics 51, 093504 (2010).
  * Mathur  _et al._ [2011] M. Mathur, I. Raychowdhury, and T. P. Sreeraj, Journal of Mathematical Physics 52, 113505 (2011).
  * Raychowdhury [2013] I. Raychowdhury,  _Prepotential Formulation of Lattice Gauge Theories_ , Ph.D. thesis, Calcutta U. (2013).
  * Zohar and Burrello [2015] E. Zohar and M. Burrello, Physical Review D 91 (2015), 10.1103/physrevd.91.054506.
  * Zohar  _et al._ [2017] E. Zohar, A. Farace, B. Reznik, and J. I. Cirac, Physical Review A 95 (2017), 10.1103/physreva.95.023604.
  * Alexandru  _et al._ [2019] A. Alexandru, P. F. Bedaque, S. Harmalkar, H. Lamm, S. Lawrence, and N. C. Warrington (NuQS), Phys. Rev. D 100, 114501 (2019), arXiv:1906.11213 [hep-lat] .
  * Ji  _et al._ [2020] Y. Ji, H. Lamm, and S. Zhu, Physical Review D 102 (2020), 10.1103/physrevd.102.114513.
  * Hackett  _et al._ [2019] D. C. Hackett, K. Howe, C. Hughes, W. Jay, E. T. Neil, and J. N. Simone, Physical Review A 99 (2019), 10.1103/physreva.99.062341.
  * Kreshchuk  _et al._ [2020] M. Kreshchuk, W. M. Kirby, G. Goldstein, H. Beauchemin, and P. J. Love, “Quantum simulation of quantum field theory in the light-front formulation,” (2020), arXiv:2002.04016 [quant-ph] .
  * Kreshchuk  _et al._ [2021] M. Kreshchuk, S. Jia, W. M. Kirby, G. Goldstein, J. P. Vary, and P. J. Love, Entropy 23, 597 (2021).
  * Buser  _et al._ [2021] A. J. Buser, H. Gharibyan, M. Hanada, M. Honda, and J. Liu, JHEP 09, 034 (2021), arXiv:2011.06576 [hep-th] .
  * Kaplan  _et al._ [2003] D. B. Kaplan, E. Katz, and M. Unsal, JHEP 05, 037 (2003), arXiv:hep-lat/0206019 .
  * Carena  _et al._ [2021] M. Carena, H. Lamm, Y.-Y. Li, and W. Liu, Phys. Rev. D 104, 094519 (2021).
  * Clemente  _et al._ [2022] G. Clemente, A. Crippa, and K. Jansen, (2022), arXiv:2206.12454 [hep-lat] .
  * Carena  _et al._ [2022] M. Carena, E. J. Gustafson, H. Lamm, Y.-Y. Li, and W. Liu, (2022), arXiv:2208.10417 [hep-lat] .
  * Preskill [2018] J. Preskill, Quantum 2, 79 (2018).