# Strategies for the Determination of the Running Coupling of
(2+1)-dimensional QED with Quantum Computing

Giuseppe Clemente<EMAIL_ADDRESS>Deutsches Elektronen-Synchrotron
(DESY), Platanenallee 6, 15738 Zeuthen, Germany  Arianna Crippa
<EMAIL_ADDRESS>Deutsches Elektronen-Synchrotron (DESY), Platanenallee
6, 15738 Zeuthen, Germany  Institüt für Physik, Humboldt-Universität zu
Berlin, Newtonstr. 15, 12489 Berlin, Germany  Karl Jansen<EMAIL_ADDRESS>Deutsches Elektronen-Synchrotron (DESY), Platanenallee 6, 15738 Zeuthen,
Germany

###### Abstract

We propose to utilize NISQ-era quantum devices to compute short distance
quantities in $(2+1)$-dimensional QED and to combine them with large volume
Monte Carlo simulations and perturbation theory. On the quantum computing
side, we perform a calculation of the mass gap in the small and intermediate
regime, demonstrating, in the latter case, that it can be resolved reliably.
The so obtained mass gap can be used to match corresponding results from Monte
Carlo simulations, which can be used eventually to set the physical scale. In
this paper we provide the setup for the quantum computation and show results
for the mass gap and the plaquette expectation value. In addition, we discuss
some ideas that can be applied to the computation of the running coupling.
Since the theory is asymptotically free, it would serve as a training ground
for future studies of QCD in $(3+1)$-dimensions on quantum computers.

##  I Introduction

In the last decades, non-perturbative numerical investigations of quantum
field theories using Markov Chain Monte Carlo (MCMC) techniques have reached
unprecedented levels of reliability and accuracy in the determination of
quantities of physical interest [1]. On the other hand, the recent advances on
quantum hardware open up the possibility of exploring new techniques and using
the Hamiltonian formulation have the potential to address new problems such as
chemical potentials, topological terms or real time evolution. However, there
are still many challenges that one has to face for both approaches. In MCMC
simulations, autocorrelation times grow rapidly, in some cases even
exponentially, (a phenomenon called critical slowing down) as the lattice
spacing vanishes, making it hard to investigate the continuum limit.
Furthermore, the path integral used for simulating some models is afflicted by
the infamous sign problem [2, 3]. In lattice quantum chromodynamics (LQCD),
this prevents, e.g., to accurately characterize the region of the phase
diagram at finite baryon density, which is of phenomenological interest [4, 5,
6].

Regarding quantum computation, with the present hardware technology only very
small systems can be simulated. However, no sign problem is present in the
Hamiltonian formalism, and no apparent obstacle hinders the investigation of
finer lattice spacings thus avoiding large autocorrelation times. For these
reasons, we believe that in the current noisy intermediate-scale quantum
(NISQ) era [7] we can use quantum computing to study small scale properties of
lattice field theories.

The system considered in this work is quantum electrodynamics (QED) in 2+1
dimensions, which has been investigated in literature with different
techniques [8, 9, 10, 11, 12, 13]. Our interest is in the study of small
distance quantities, a most prominent example being the running coupling. By
making contact with perturbation theory this allows to calculate also the
$\Lambda$-parameter, $\Lambda^{QED}$ in 2+1-dimensional QED, analogous to the
QCD $\Lambda$-parameter. Other possible quantities of interest that is
possible to study are, for example, the (scale dependent) renormalization
factors $Z$ [14, 15].

It is the main idea of this and follow-up works, to match the results of short
distance quantities obtained from quantum computations with the ones coming
from MCMC simulations in the strong and intermediate coupling region. In
particular, when such a matching is successfully carried out, large scale MCMC
calculations will provide us with a physical value of the lattice spacing.
This, in turn, allows us to convert results which are obtained in lattice
units, e.g. the renormalization scale or the $\Lambda$-parameter, to physical
units.

The paper is organized as follows. In Sec. II we introduce the system under
investigation: the lattice discretization and Hamiltonian formulation are
described in Sec. II.1, while the electric and magnetic bases are discussed in
Sec. II.2. In Sec. III we introduce the numerical setup, in particular, the
encoding adopted (Sec. III.1) and the variational technique (Sec. III.2).
Numerical results are discussed in Sec. IV, where we show the behavior of the
spectral gap and plaquette observable with both exact diagonalization data and
the results of the application of variational quantum techniques on a
simulator. In Sec. V we formulate some proposal to match quantum results with
Monte Carlo simulations and to compute the running of the coupling in a non-
perturbative way. Finally, in Sec. VI we summarize our findings, and discuss
future perspectives.

Two proposal for the improvement of the ansatz and for a more efficient
encoding are considered in Appendix A and B respectively, while Appendix C is
devoted to present tests and results for the penalty term method performed
with the pure gauge system and the fermionic case.

##  II (2+1)-dimensional QED on the lattice

As for QCD in 3+1 dimensions, the behavior of QED in 2+1 dimensions presents
confinement and asymptotic freedom. In this section we describe the
discretization adopted [15].

###  II.1 Hamiltonian

In order to deal with the fermionic doubling problem [16, 17, 18], i.e. the
existence in $d$-dimensions of $2^{d}$ flavors for each physical particle, we
consider a Kogut-Susskind (K-S) formulation [19].

$(0,0)$$(0,1)$$(1,1)$$(1,0)$$\hat{E}_{\vec{n},y},\hat{U}_{\vec{n},y}$$\hat{E}_{\vec{n},x},\hat{U}_{\vec{n},x}$
Figure 1: One plaquette system. Fermions are on the $\vec{n}$ site, fields and
the link operators $\hat{U}$ on the links. The arrows follow the positive
directions $x,y$.

With this formulation the fermionic degrees of freedom are discretized on
distinct lattices and then separated by $2a$. Fermions and antifermions
formulation are represented by single component field operators
$\hat{\phi}_{\vec{n}}$, with $\vec{n}=(n_{x},n_{y})$ as the position on the
lattice sites, as in Figure 1.

The Hamiltonian that we use to represent QED [20, 21] consists of four terms:

$H_{tot}=H_{E}+H_{B}+H_{m}+H_{kin}$ (1)

The first term is

$\displaystyle\hat{H}_{E}=\frac{g^{2}}{2}\sum_{\vec{n}}\left(\hat{E}^{2}_{\vec{n},x}+\hat{E}^{2}_{\vec{n},y}\right),$
(2)

where the field operators $\hat{E}_{\vec{n},\mu}$ (direction $\mu=x,y$) are
dimensionless and act on the links, and $g/\sqrt{a}\to g$ is the coupling
constant. The electric fields take integer eigenvalues $e_{\vec{n}}=0,\pm
1,\pm 2,...$ with
$\hat{E}_{\vec{n},{\mu}}\ket{e_{\vec{n}}}=e_{\vec{n}}\ket{e_{\vec{n}}}$.

In LGTs, the following Wilson operators are introduced [22] on the links,

$\displaystyle\hat{U}=\mathrm{e}^{iag\hat{A}_{\mu}(\vec{n})},$ (3)

where $\hat{A}_{\mu}(\vec{n})$ is the vector field. They measure the phase
proportional to $g$ acquired by a unit charge moved along the link and act as
a lowering operator on electric field eigenstates:
$\hat{U}_{\vec{n},{\mu}}\ket{e_{\vec{n}}}=\ket{e_{\vec{n}}-1}$. Thus, we have
the following commutator

$[\hat{E}_{\vec{n},\mu},\hat{U}_{\vec{n}^{\prime},\nu}]=-\delta_{\vec{n},\vec{n}^{\prime}}\delta_{\mu,\nu}\hat{U}_{\vec{n},\mu}.$
(4)

The magnetic interaction is determined by the plaquette operator
$\hat{P}_{\vec{n}}=\hat{U}_{\vec{n},x}\hat{U}_{\vec{n}+x,y}\hat{U}^{{\dagger}}_{\vec{n}+y,x}\hat{U}^{{\dagger}}_{\vec{n},y}$,

$\displaystyle\hat{H}_{B}=-\frac{1}{2g^{2}}\sum_{\vec{n}}\left(\hat{P}_{\vec{n}}+\hat{P}_{\vec{n}}^{{\dagger}}\right).$
(5)

For the fermionic part of the Hamiltonian we have the mass term

$\displaystyle\hat{H}_{m}=m\sum_{\vec{n}}(-1)^{n_{x}+n_{y}}\hat{\phi}^{\dagger}_{\vec{n}}\hat{\phi}_{\vec{n}},$
(6)

and the kinetic term,

$\displaystyle\hat{H}_{kin}=\Omega\sum_{\vec{n}}\sum_{\mu=x,y}(\hat{\phi}_{\vec{n}}^{{\dagger}}\hat{U}_{\vec{n},\mu}\hat{\phi}_{\vec{n}+\mu}+{H.c.}),$
(7)

where $[\Omega]=[1/(2a)]$. It corresponds to the creation or annihilation of a
fermion-antifermion pair on neighbouring lattice sites and the adjustment of
the link.

The Hamiltonian is gauge invariant, i.e. it commutes with the Gauss’ law
operators $G_{\vec{n}}$ at each site $\vec{n}$

$\displaystyle G_{\vec{n}}=\Bigg{[}\sum\limits_{\mu=x,y}$
$\displaystyle\left(\hat{E}_{\vec{n},\mu}-\hat{E}_{\vec{n}-\mu,\mu}\right)-\hat{q}_{\vec{n}}-\hat{Q}_{\vec{n}}\Bigg{]}$
(8) $\displaystyle
G_{\vec{n}}\ket{\Phi}=0\iff\ket{\Phi}\in\mathcal{H}_{\text{phys.}},$

where $\hat{q}_{\vec{n}}$ are dynamical charges, defined as

$q_{\vec{n}}=\hat{\phi}^{\dagger}_{\vec{n}}\hat{\phi}_{\vec{n}}-\frac{1}{2}(1+(-1)^{\lvert\vec{n}\rvert+1}).$
(9)

$\hat{Q}_{\vec{n}}$ are external static charges, which are set to non-zero
values only when one considers the computation of the static potential.

###  II.2 Electric and magnetic basis and fermion discretization

In this section we consider the calculation of the energy gap with periodic
boundary conditions (PBC). These allow, unlike open boundary conditions, to
properly define the momentum operator and hence project into the zero-momentum
sector ($\vec{p}_{\text{tot.}}=\vec{0}$) along with the total zero-charge
sector ($Q_{\text{tot.}}=0$). In principle it should be easier to extract the
energy difference when the number of states to be considered is smaller.
However, this procedure is not straightforward in the open boundary case,
since no momentum operator in the sense of generator of translations exists;
in this case, it is still possible to introduce a pseudo-momentum operator
[23] whose square can be used to select only the physical states among all the
results a posteriori, or as a suppression term with the same procedure
described in Sec. III.3.

$R1$$R2$$R3$$R4$$R_{y}$$R_{x}$ Figure 2: $2\times 2$ PBC lattice. By applying
the Gauss’ law we can write $R4$ (thin line) in terms of the other three
rotators.

We consider a square lattice formed by four lattice sites. Due to Gauss’ law,
only five of the eight links are independent, thus we can have a more resource
efficient Hamiltonian, suitable for a wide range of quantum hardware. To
conveniently simplify the expressions we use the set of operators, i.e.
rotators and strings [20], which preserve Gauss’ law, as in Figure 2. Using
periodic boundary conditions, the momentum-zero sector can be enforced exactly
by using a translational invariant Hamiltonian111 However, after gauge fixing
and for a finite truncation, symmetry cannot be enforced exactly on the
remaining gauge degrees of freedom (i.e., rotators and strings)..

As mentioned, our analysis include fermionic matter and excited states and by
using a K-S formulation we have a symmetric fermionic system every two lattice
spacings.

####  II.2.1 Discretization and Truncation

In the following, we describe two schemes that allow to discretize the
continuous $U(1)$ group with $\mathbb{Z}_{2L+1}$ which provides a discrete
basis for the operators $\hat{A}_{\vec{n},{\mu}}$. Since the gauge rotators
possess discrete but infinite spectra, any numerical approach requires a
truncation of the Hilbert space. Thus the gauge fields assume integer values
in range $[-l,l]$, and the dimension of the Hilbert space (within the gauge
part) is $(2l+1)^{5}$. The parameter $l$ is the cut-off value for the group
truncation.

We consider two different representations for the Hamiltonian, i.e. electric
and magnetic. While in the former the operators associated to the electric
field are diagonal, the plaquette operator corresponds to a non-diagonal
expression. As a consequence, applying a different representation in the weak
coupling regime is preferable. This is performed through the discrete Fourier
transform, that diagonalizes the lowering operators.

The discrete approximation in the magnetic basis is related to the parameters
$l$ and $L$. We refer to $l$ as the truncation level and $L$ as the
discretization level. Since the truncated $U(1)$ and $\mathbb{Z}_{2L+1}$ are
equivalent in the electric basis, the parameter $L$ is irrelevant, however it
strongly influences the results derived in the magnetic representation. For a
given $l$, $L$ defines the resolution of the approximation centered around the
vacuum. The group is discretized into $2L+1$ states and only $2l+1$ are kept
after the truncation.

####  II.2.2 Jordan-Wigner

For the numerical implementation, the fermionic degrees of freedom can be
mapped to spins using a Jordan-Wigner transformation [24],

$\displaystyle\phi_{\vec{n}}$
$\displaystyle=\Big{[}\prod\limits_{{\vec{k}}<{\vec{n}}}(i\sigma^{z}_{\vec{k}})\Big{]}\sigma^{-}_{\vec{n}},$
(10) $\displaystyle\phi_{\vec{n}}^{\dagger}$
$\displaystyle=\Big{[}\prod\limits_{{\vec{k}}<{\vec{n}}}(-i\sigma^{z}_{\vec{k}})\Big{]}\sigma^{+}_{\vec{n}},$
(11)

where $\sigma^{z}$ is the $z-$Pauli matrix,
$\sigma^{\pm}\equiv\frac{\sigma^{x}\pm i\sigma^{y}}{2}$ and the relation
$\vec{k}<\vec{n}$ is defined by $(0,0)<(0,1)<(1,1)<(1,0)$ to satisfy the
fermionic commutation relations.

Applying Eq. (10) to (9), the expressions for the dynamical charges in terms
of bosonic spin operators become:

$q_{\vec{n}}=\frac{1}{2}(\sigma_{\vec{n}}^{z}+(-1)^{\lvert\vec{n}\rvert+1}).$
(12)

##  III Numerical setups and methods

###  III.1 Encoding

The encoding adopted in [21] maps the $N$ fermionic states into an equal
number of qubits and gauge physical states onto $2l+1$ qubits using

$\displaystyle\ket{-l+j}_{\text{phys.}}$
$\displaystyle\mapsto\ket*{\overbrace{0\ldots 0}^{j}1\overbrace{0\ldots
0}^{2l-j}}.$ (13)

The truncated electric field and link operators (for a single edge) can be
expressed starting from their expressions in the physical basis:

$\displaystyle\hat{E}$
$\displaystyle=\sum_{i=-l}^{l}i\ket{i}_{\text{phys.}}\bra{i}_{\text{phys.}},$
(14) $\displaystyle\hat{U}$
$\displaystyle=\sum_{i=-l+1}^{l}\ket{i-1}_{\text{phys.}}\bra{i}_{\text{phys.}}.$
(15)

For example, for $l=1$ the corresponding $2l+1$ possible states are encoded as

$\displaystyle\ket{-1}_{\text{phys.}}\mapsto$ $\displaystyle~{}\ket{100},$
(16) $\displaystyle\ket{0}_{\text{phys.}}\mapsto$ $\displaystyle~{}\ket{010},$
(17) $\displaystyle\ket{1}_{\text{phys.}}\mapsto$ $\displaystyle~{}\ket{001}.$
(18)

However, this mapping is not resource efficient since we would need $2l+1$
qubits for each gauge variable.

Hereby we use another method to encode the physical states. Instead of using
$2l+1$ qubits for each truncated gauge variable, one would like to use a
compact encoding, see e.g. [25]. The minimum number of qubits required per
gauge variable is222In principle, one could even encode in the qubit space any
number $N$ of gauge variables, so that the space one uses
$\lceil\log_{2}[N(2l+1)]\rceil$ qubits instead of
$N\lceil\log_{2}(2l+1)\rceil$, gaining some additional space, but not much.
$q_{\text{min}}=\lceil\log_{2}(2l+1)\rceil$, such that, using $N$ truncated
gauge variables, the fraction of physical states ${(2l+1)}^{N}$ among the
total number of states in the encoding is $2^{N(\\{\log_{2}(2l+1)\\}-1)}$,
which corresponds to $\frac{1}{2^{N}}{(1-\frac{1}{4l})}^{N}$ in the worst case
(when $l$ is a power of $2$); using the previous encoding instead, one would
use just a fraction $\frac{{(2l+1)}^{N}}{2^{N(2l+1)}}$, which vanishes
exponentially in the limit $l\to\infty$ with a rate coefficient $2N\ln 2$.

In order to represent a generic transition in the following discussion, we use
the notation

$\displaystyle\sigma^{+}$ $\displaystyle\equiv\outerproduct{1}{0},\,$
$\displaystyle\sigma^{\uparrow}$ $\displaystyle\equiv\outerproduct{1}{1},\,$
(19) $\displaystyle\sigma^{-}$ $\displaystyle\equiv\outerproduct{0}{1},\,$
$\displaystyle\sigma^{\downarrow}$ $\displaystyle\equiv\outerproduct{0}{0}.$

Amongst all the possible $\binom{2^{q_{\text{min}}}}{2l+1}$ compact encodings,
it is convenient to pick the ones that simplify the ladder terms
$\ket{i-1}_{\text{phys.}}\bra{i}_{\text{phys.}}$ in Eq. (15) in such a way
that the bit string representing the state in the computational basis changes
only by one bit. Another way of phrasing this condition is that we would
minimize the number of off-diagonal operators appearing in the Pauli string
decompositions of these terms. This is realized by the so called Gray codes.

For example, with $l=1$ there are three physical states
$\ket{i}_{\text{phys.}}$ for $i\in\\{-1,0,1\\}$, which can be encoded using
only $2$ qubits using a Gray encoding pattern as shown in Table 1, where the
state $\ket{10}$ is considered unphysical.

Table 1: Gray encoding for $l=1$. $\ket{i}_{\text{phys.}}$ | $\ket{i}$ | $\outerproduct{i}{i}$ | $\outerproduct{i-1}{i}$
---|---|---|---
$\ket{-1}_{\text{phys.}}$ | $\ket{00}$ | $\sigma^{\downarrow}\otimes\sigma^{\downarrow}$ | $0$
$\ket{0}_{\text{phys.}}$ | $\ket{01}$ | $\sigma^{\downarrow}\otimes\sigma^{\uparrow}$ | $\sigma^{\downarrow}\otimes\sigma^{-}$
$\ket{+1}_{\text{phys.}}$ | $\ket{11}$ | $\sigma^{\uparrow}\otimes\sigma^{\uparrow}$ | $\sigma^{-}\otimes\sigma^{\uparrow}$
unphysical | $\ket{10}$ |  |

The expressions for the truncated electric field and link operators then
become:

$\displaystyle\hat{E}$
$\displaystyle\mapsto-\outerproduct{00}{00}+\outerproduct{11}{11}=-\frac{1}{2}\Big{[}\sigma^{z}_{0}+\sigma^{z}_{1}\Big{]},$
(20) $\displaystyle\hat{U}$
$\displaystyle\mapsto\outerproduct{00}{01}+\outerproduct{01}{11}=\frac{1}{2}\Big{[}\sigma^{-}_{0}(I_{1}+\sigma^{z}_{1})+\sigma^{-}_{1}(I_{0}-\sigma^{z}_{0})\Big{]}.$
(21)

Unfortunately, since $2l+1$, the number physical basis states for gauge
variables, is not a power of two, one has to deal with the unphysical states
so that they would not appear in the final solutions. This can be done either
by working with a structured ansatz with rotations allowed only between
physical states (see Appendix A) or by introducing a penalty term in the
Hamiltonian which suppress solutions with overlap in the unphysical region of
the Hilbert space (see Sec. III.3). Another approach which can be applied in
the case of magnetic basis consists in changing the discretization in such a
way that one can include the remaining $2^{\lceil\log_{2}(2l+1)\rceil}-2l-1$
states as physical (see Appendix B).

###  III.2 Variational Quantum Algorithms

In this section we describe the technique adopted to compute the mass gap and
the expectation value of the plaquette operator, using an extension of the
Variational Quantum Eigensolver (VQE) algorithm [26]. A quantum device can be
used to efficiently evaluate the expectation value of tensor products of an
arbitrary sequence of Pauli operators [27]. Since the Hamiltonian is written
as a weighted sum of such Pauli string terms, also $\langle H\rangle$ can be
efficiently estimated333In general, it is possible to choose a grouping
strategy to identify subsets of Pauli strings appearing in the Hamiltonian, it
is possible to reduce the number of independent circuit evaluations [28]..

In order to find the lowest eigenvalue of a given operator $H$, the
variational approach finds an approximation to the eigenvector $\ket{\psi}$
which corresponds to the lowest eigenvalue and that minimizes

$E(\theta):=\expectationvalue{H}{\psi(\theta)},$ (22)

where the state $\ket{\psi({\theta})}=U(\theta)\ket{0}$ is realized as a
parameterized circuit $U(\theta)$ called ansatz.

This procedure is done by varying a vector $\theta$ of scalar parameters
(typically gate rotation angles) through the combination of a classical and a
quantum part.

The Variational Quantum Deflation (VQD) method [29, 30] extends VQE to
estimate the $k$-th excited state $E_{k}$ by penalizing the solutions of the
lowest excited states. This is done through a minimization of the cost
function

$C(\theta_{k})=\expectationvalue{H}{\psi(\theta_{k})}+\sum_{i=0}^{k-1}\beta_{i}{\big{\lvert}\innerproduct{\psi(\theta_{k})}{\psi(\theta_{i}^{*})}\big{\rvert}}^{2},$
(23)

where $\beta_{i}$ are real-valued coefficients (which must be larger than the
gaps $E_{k}-E_{i}$) and $\theta_{i}^{*}$ are the optimal parameters for the
$i$-th excited state. The overlap terms are computed by either using the
inverse circuit444Here and in our results we assume the parameterized ansatz
to be the same $U(\theta)$ for every excited state, but this requirement is
not necessary. $U^{\dagger}(\theta_{k})U(\theta_{i}^{*})\ket{0}$ and
estimating the occurrence of all-zero measurements or by a SWAP test [30, 31,
32]. This can be interpreted as minimizing the $E(\theta_{k})$ with the
constraint that $\ket{\psi(\theta_{k})}$ must be orthogonal to the previous
$k$ states. Since our goal is to compute the energy gap between the ground
state $E_{0}$ and first excited state $E_{1}$, we follow three main steps:

  * 1.

Perform the VQE and obtain optimal parameters and an approximate ground state
$\ket{\psi(\theta_{0}^{*})}$;

  * 2.

For $E_{1}$ define a Hamiltonian:

$H_{1}=H+\beta\ket{\psi(\theta_{0}^{*})}\bra{\psi(\theta_{0}^{*})};$ (24)

  * 3.

Perform the VQE with the Hamiltonian $H_{1}$ to find an approximation of the
first excited state $\bra{\psi(\theta_{1}^{*})}$.

In our case, for each value of the coupling $g$ we can compute the best
approximation to the ground state and extract $E_{0}$ and the expectation
value of the plaquette $\expectationvalue{\Box}$, while from the first excited
state we just need $E_{1}$ to estimate the spectral gap $\Delta
E=E_{1}-E_{0}$.

###  III.3 Ansatz and penalty term

In this section we describe the procedure that we used to deal with unphysical
states. Instead of constraining the reachable states to the physical ones at
the level of the ansatz (as discussed in Appendix A for the gray encoding), we
consider a generic ansatz and introduce, directly in the definition of the
Hamiltonian, a penalty term that suppresses unphysical contributions on the
final states [10]. The form that we used for this suppression term is the
following:

$\Delta
H_{\text{suppr.}}=\lambda\Big{[}\sum\limits_{p\in{\text{gauge}}}\Pi^{(u)}_{p}+\Pi^{(Q_{\text{tot}}\neq
0)}_{\text{ferm.}}\Big{]},$

where $\lambda$ is the suppression coefficient, while $\Pi^{(u)}_{p}$ and
$\Pi^{(Q_{\text{tot}}\neq 0)}_{\text{ferm.}}$ are respectively the projectors
onto the unphysical Hilbert spaces of single gauge variables (identity on the
other variables) and the projector onto the non-zero charge space.

At the end of the VQE optimization, we can assess how much the optimal state
reached $\ket{\psi(\theta^{*})}$ is unphysical by measuring its overlap with
the unphysical Hilbert space $\mathcal{H}_{\text{unphys.}}$. In practice, this
is done by computing the expectation value of the projector into
$\mathcal{H}_{\text{unphys.}}$,

$u(\theta^{*}_{k})\equiv\expectationvalue{\Pi_{\text{unphys.}}}{\psi(\theta^{*}_{k})}.$
(25)

First, we tested this method for the pure gauge system, then we considered the
fermionic case with a selection procedure for $\lambda$ based on a bisection
approach. See Appendix C for details and results.

##  IV Numerical results

Here we discuss results of exact diagonalization (ED) of the Hamiltonian Eq.
(1) in both electric and magnetic basis and with $\Omega=1$, $m=0$. Then we
show some results of the VQD method.

###  IV.1 Exact diagonalization results

Figure 3 and Figure 4 show respectively the ground state energy and plaquette
expectation value for different levels of truncation $l$, and for the
discretization level $L=8$ in the magnetic basis. As expected, for small
values of the coupling $g$ the convergence in the truncation parameter $l$ at
fixed discretization $L$ is generally faster in the magnetic basis. In the
case of larger couplings the electric basis performs better, due to the
relative dominance of the $H_{E}$ or $H_{B}$ terms in the Hamiltonian.

Figure 3: Exact diagonalization data for ground state energy in both electric
and magnetic basis and selected values of discretization and truncation (see
legend). Figure 4: Exact diagonalization data for the plaquette in both
electric and magnetic basis and selected values of discretization and
truncation (see legend). Figure 5: Exact diagonalization data for ground state
energy as a function of $1/l$ for both electric and magnetic basis and
selected values of $g$ (see legend). Figure 6: Exact diagonalization data for
spectral gap in both electric and magnetic basis and selected values of
discretization and truncation (see legend).

This behavior is apparent in Figure 5, where we plot the estimates to the
untruncated theory as functions of the inverse truncation $1/l$ for some
selected values of the coupling and for two values of $L$, i.e. 4 and 8. For
the electric basis, the limit $l\to\infty$ is sufficient to discuss about
convergence to the untruncated theory, while for the magnetic basis one has to
consider a double limit where also the discretization level $L$ is taken into
account. However, analyzing the spectral gap in the range $g\in[0.2,0.8]$,
shown in Figure 6, one sees a poorer convergence in $l$ of the magnetic basis
results with the electric ones, which for $l=3$ are already quite close to the
best estimate of the untruncated theory, which is illustrated by the dashed
line and is represented by the electric basis result at $l=9$. For smaller
values of $g$ we expect the magnetic basis to be more convenient, but since
the matching with MCMC results can be done for the couplings of $O(1)$ or
larger, we decided to consider mainly the electric basis for the VQD results
in the following discussion.

###  IV.2 VQD results

For the VQD results we adopted the NFT optimizer [33] with the Qiskit’s
EfficientSU(2) [34] generic ansatz up to $5$ layers and with full entanglement
gates (i.e., $CNOT$ gates between every pair of qubits) alternated to single-
qubit rotation layers. Since we are testing the feasibility of the method,
numerical results in this paper have been obtained using a simulator without
noise and with an infinite number of shots.

We first launched some runs to determine the best estimate to the ground state
and its energy. The results of the ground state energy and the plaquette as a
function of $g^{-2}$ in the range $[0.2,3]$ for $l=1,2,3$ in the electric
basis are shown in Figure 7 and 8. Then we used the best ground state
obtained, for each $g$ and $l$, as a further penalty term in order to find the
energy of the first excited state, as described in Sec. III.2. From the exact
result we have seen that the gap is quite small, thus we must be careful
during the selection of the $\lambda$ factor and the number of iterations. As
mentioned in Sec. III.3,we consider a bisection method in order to choose the
best value of $\lambda$. Starting from an initial value of $\lambda\sim
O(1)$555The factor must be at least larger than the energy gap, as depicted in
the pure gauge plots in Appendix C.1., we compute the percentage of unphysical
state and tune the suppression factor in terms of this percentage. In
particular we choose a threshold of $99\%$ of physical component.

Figure 7: Best results for VQD ground state energy as a function of the
coupling in the electric basis (dots) at some values of truncation level $l$
and exact diagonalization (lines). Bottom panel: discrepancies with the exact
values.

Regarding the expectation value of the plaquette, best results of the VQE
(which minimize $E_{0}$) are shown in Figure 4. The discrepancy between VQE
and exact diagonalization results for $\expectationvalue{\Box}$ is larger
especially in the region close to $g=1$. This may be explained by two
competing factors. On one hand, due to the closing of the gap (see Figure 6)
and the accumulation of quasi-degenerate levels in the spectrum for smaller
values of $g$, the quality of the ground state obtained
$\ket{\psi(\theta_{0}^{*})}$ (and therefore also the estimate of the
plaquette) is affected, since it is more likely to converge to a superposition
of quasi-degenerate states which cannot be efficiently discriminated by the
$VQE$ optimization process. On the other hand, besides an irrelevant factor,
the plaquette observable coincides with the magnetic part of the Hamiltonian,
which becomes dominant in the regime of small $g$, so that a decent estimate
of the ground state energy $E_{0}$ (which is exactly what the VQE optimizes)
coincides with an acceptable estimate also for $\expectationvalue{\Box}$.

Figure 8: Plaquette measurements on the ground state in Figure 7 as a function
of the coupling $g$ ($m=0$) in the electric basis. VQE results (dots) and ED
(lines).

The ED and best VQD results of the spectral gap computed in the electric basis
are shown in Figure 9 for truncation level up to $l=3$. As mentioned
previously, the closing of the gap for $g<1$, corresponding to the shaded area
on the right of Figure 9, makes it difficult for the VQD to reach enough
accuracy in the estimation of both the ground state and first excited levels.
In the region $g\geq 1$ the convergence is more under control. We conclude
that a matching between MCMC data and VQD using the mass gap as observable
could be possible only in the region $g\geq 1$, because, due to physical
features of the lower part of the spectrum, the accuracy required does not
allow to realistically extend the matching to smaller values of $g$. However,
as we describe in the next section, in order to extend the study of the
running coupling to the weak coupling regime with $g<1$, we would use
information from the static force instead of the mass gap.

Figure 9: Best results for spectral gap as a function of the coupling $g$
($m=0$) in the electric basis. VQD results (dots) and ED (lines); the shaded
area corresponds to the region where we could not obtain results with enough
precision to estimate the gap reliably. For all data points shown, the
unphysical part of both the ground state and first excited does not exceed
$0.5\%$.

Since the ground state and first excited levels are well separated at values
of the coupling larger than $g=1$, the first VQE optimization stage succeeds
with higher probability in finding a good approximation
$\ket{\psi(\theta^{*}_{0})}$ to the ground state. For values of the coupling
smaller than $g=1$ the ground state and excited levels get closer and closer,
making the optimization process harder, since these states become quasi-
degenerate and more iterations are required to distinguish the optimal
direction from the others in the lower spectrum. This can also be seen in the
bottom panel of Figure 7, which shows the discrepancy between the VQE results
and the exact values.

##  V Outlook: towards the running coupling

As outlined in the previous sections, it is the main purpose of this work to
establish the connection between quantum computations of short distance
quantities and large volume Monte Carlo simulations which can provide us with
the value of the lattice spacing. To this end, we have described here the
setup for the quantum computation and showed that we can compute the mass gap
in the relevant region of bare couplings where we expect to be able to match
Monte Carlo simulations.

As an important example of a short distance quantity we consider the running
coupling. Despite the limited computational resources of the present quantum
computer hardware, it is still possible to provide definitions of the running
coupling using the plaquette observable, the static potential or the static
force at short distances. In the following, we describe both approaches
shortly, devoting a more detailed discussion together with first numerical
results to a future work.

###  V.1 Matching with MCMC data

To set the scale, i.e., the physical values of the lattice spacing, we propose
to take into consideration the spectral gap $\Delta E$ in the scalar sector
and the static potential $V(r)$, which is defined as the lowest energy of a
system of a static quark-antiquark pair at distance $r$.

As in QCD, from both lattice and experimental evaluations of $V(r)$, up to an
unphysical offset, one can determine the analogue of the _Sommer parameter_
$r_{0}$ using the force [35]:

$r_{0}^{2}\frac{\partial V(r)}{\partial r}\Big{\rvert}_{r=r_{0}}=c,$ (26)

where $c$ is an suitably chosen constant which in QCD is conventionally fixed
to $c^{QCD}=1.65$ with a phenomenological value $r_{0}^{QCD}=0.5$ fm. In the
continuum, the expected form of the static potential in 2+1-dimensional QED is
the following [36, 37, 38]

$V(r)=V_{0}+\alpha\log r+\sigma r,$ (27)

where the second term is a logarithmic Coulomb term while the third one
represents the confinement term with string tension $\sigma$. On the lattice,
can discretize the derivative in Eq. (26) by finite differences
$\frac{\partial V(r)}{\partial r}\simeq\frac{V(r_{2})-V(r_{1})}{r_{2}-r_{1}}$.
Despite the current size limitations for lattices on quantum machines, it is
sufficient to evaluate $V(r)$ using just $2$ distances, which means $3$
lattice locations where the static charges are placed. As we describe in more
details in the following section, in general the QED results for a single
fermionic species depend on both the bare gauge coupling $g$ and the bare
fermionic mass $m$. From the static potential computed on a lattice one can
estimate the static force at two distances $r$ and as a function of the
couplings $F(r;g,m)$. The force and the mass gap can be computed both, from
MCMC simulations and in the Hamiltonian formalism. In order to match both
approaches, the task is then to find matching values $(g^{*},m^{*})$ where the
results for the mass gap the force agree. Such a matching is expected to be
possible in the intermediate coupling regime where, as we show in this paper
with our Hamiltonian setup, we can reach sufficiently accurate results for the
mass gap. A numerical investigation and matching of the static force will be
discussed in a future work.

###  V.2 Step scaling approach

One approach for computing the running coupling is to “scale step” [39, 40],
see [41] for an introduction to the Schrödinger functional scheme for the step
scaling procedure.

As an example, we consider the coupling defined from the static potential. Let
us take, for illustration purposes a Coulomb potential666The here discusses
example can be straightforwardly extended to the case of 2+1-dimensional QED
with the potential given by Eq. (27). $V(r)=\alpha(r)/r$ and a pure gauge
theory with bare coupling $g$. Then, a running coupling at two scales $r_{1}$
and $r_{2}$ can be defined by $\alpha_{\mathrm{r}en}(r_{1})=r_{1}V(r_{1})$ and
$\alpha_{\mathrm{r}en}(r_{2})=r_{2}V(r_{2})$. Let us now assume that we have
computed $\alpha_{\mathrm{r}en}(r_{1}^{i=0},g_{0})$ and
$\alpha_{\mathrm{r}en}(r_{2}^{i=0},g_{0})$ at a fixed value of the bare
coupling $g^{2}=g_{0}^{2}$ at a step $i=0$. The two distances can be related
by a scale factor $s$, e.g. $s=2$ such that $r_{2}^{i=0}=sr_{1}^{i=0}$. In the
next step, $i=1$, the bare coupling $g^{2}$ is tuned in such a way that
$\alpha_{\mathrm{r}en}(r_{1}^{i=1},g)=\alpha_{\mathrm{r}en}(r_{2}^{i=0},g_{0})$
which would provide a value $g_{1}^{2}$ where the renormalized couplings agree
and hence the scales match. Applying the same scale factor $s$ at the found
bare coupling $g_{1}^{2}$ one arrives at the renormalized coupling at the
scale $r_{2}^{i=1}=sr_{1}^{i=1}=2sr_{1}^{i=0}$. Thus we get sequence of
renormalized couplings $\alpha_{\mathrm{r}en}(r_{1}^{i=0},g_{0})$,
$\alpha_{\mathrm{r}en}(2sr_{2}^{i=1},g_{1})$.

This procedure can be repeated $N$ times such that we obtain the scale
dependence of the coupling where the scale changes by the factor $s$ in each
step arriving thus at a renormalized coupling
$\alpha_{\mathrm{r}en}(Nsr_{1}^{i=0},g_{\mathrm{match}})$ with
$g_{\mathrm{match}}$ chosen large enough that contact with Monte Carlo
simulations can be made. Note that up to this point we have only worked with
quantities in lattice units thus not knowing the physical values of the scale.

It is exactly at this point where the matching of the quantum and the
classical Monte Carlo computations comes into play. The Monte Carlo simulation
can be performed in large volumes at intermediate values of the coupling and
through the strategy described below a value of the lattice spacing can be
determined. In this way, the final scale can be converted to physical units,
i.e. $r_{\mathrm{phys}}=aNsr_{1}^{i=0}$, given thus the renormalized coupling
$\alpha_{\mathrm{r}en}(r_{\mathrm{phys}},g_{\mathrm{match}})$. The sequence of
couplings obtained as outlined above can now be inverted by subsequently
changing the scale by a factor $s$, i.e. $r_{\mathrm{phys}}^{i=N}\rightarrow
r_{\mathrm{phys}}^{i=N-1}=r_{\mathrm{phys}}^{i=N}/s$ $\rightarrow
r_{\mathrm{phys}}^{i=N-2}=r_{\mathrm{phys}}^{i=N-1}/s$, …, $\rightarrow
r_{\mathrm{phys}}^{i=0}$ with the corresponding changes of the coupling. In
this way, the renormalized coupling is obtained as a function of the physical
scale and, by making contact with perturbation theory, eventually also the
important $\Lambda$ parameter, which provides the scale where non-perturbative
physics sets in, can be determined in physical units.

The procedure explained here can, in principle, also be used to disentangle
lattice effects from the real running of the coupling by taking the continuum
limit. However, this would require large lattices and, as mentioned already,
with present quantum computing hardware resources this is not feasible and
would require future quantum computers with more and improved –ideally error
corrected– qubits. Nevertheless, the just described procedure can be
implemented on already existing quantum devices allowing thus to go to the
deep perturbative regime and to make contact to low order perturbation theory.

###  V.3 Boosted coupling approach and scale setting

An alternative way to determine the running coupling is to employ the boosted
coupling defined by

$g^{2}_{\rm boosted}=g^{2}/\langle\Box\rangle\;,$ (28)

where $\expectationvalue{\Box}$ denotes the expectation value of the plaquette
operator and $g^{2}$ is the bare gauge coupling. Following the strategy given
in [14] the scale dependence of the coupling can then be determined employing
perturbation theory to a given order. An essential element in this procedure
is the determination of the renormalization scale in physical units for which
the lattice spacing needs to be calculated. The general setup for such a
calculation is illustrated in Fig. 10 where $\expectationvalue{O_{1}}$ and
$\expectationvalue{O_{2}}$ denote expectation values of two observables777The
ratio $R$ ought be dimensionless, meaning that $O_{1,\mathrm{phys}}$ and
$O_{2,\mathrm{phys}}$ have the same mass dimension. In Fig. 10 we assume for
simplicity a mass dimension of one in order to relate the lattice value
$\langle O_{1}^{\mathrm{latt}}\rangle$ to the physical one, $\langle
O_{1}^{\mathrm{phys}}\rangle$, and hence extract the lattice spacing. which
can, in principle, be extracted from experiments. Thus the physical value of
the ratio $R=\expectationvalue{O_{1}}/\expectationvalue{O_{2}}$ is known
(indicated by a ”$\bullet$” in Fig. 10). One can then tune the parameters of
the theory such that at a certain value of
$\expectationvalue{O_{1,\mathrm{latt}}}$, $R$ is reproduced. At this value of
$\expectationvalue{O_{1,\mathrm{latt}}}$ the lattice spacing can be determined
by the relation $aO_{1,\mathrm{phys}}=O_{1,\mathrm{latt}}$.

Figure 10: Illustration of principle way to determine a value of the lattice
spacing, see discussion in main text.

Examples of $O_{1}$ and $O_{2}$ are particle masses (mass gaps) or the static
force at a given physical distance. However, many more choices are possible
and actually used in large scale lattice simulations. As we show in this
paper, mass gaps become increasingly difficult to determine when the coupling
is decreased towards the continuum limit and hence we, unfortunately, consider
this not to be an option when we want to work in the perturbative regime.

As an alternative we consider the static force at small distances. In this
case, we can define $\langle O_{j,\mathrm{phys}}\rangle$ as
$F_{\mathrm{phys}}(r_{j,\mathrm{phys}})$. A difficulty for using this strategy
is that on the lattice the physical distance $r_{\mathrm{phys}}=aN$, with $N$
the number of lattice points, needs to be kept constant. This can be achieved,
e.g., by demanding that $r^{2}F(r)=\mathrm{constant}$. As a consequence, when
going to smaller values of the lattice spacing, the number of lattice points
needs to be increased correspondingly. We therefore consider the possibility
of using the static force for setting the scale as a conceptually clean way
which can, however, only be employed when significantly larger future quantum
computing resources are available. Hence, we consider the here described
strategy to compute the running coupling as a nice but only future
perspective.

##  VI Conclusions

In this work we have provided a first step to combine classical Markov Chain
Monte Carlo simulations with quantum computations. While MCMC calculations can
be performed on very large lattices nowadays, they are limited to reach small
values of the lattice spacing. On the other hand, quantum computations which
uses the Hamiltonian formulation can, at least in principle, be used at
arbitrary lattice spacings. However, presently, quantum computations are
limited to small systems given the present generation of noisy quantum
computers. They are therefore restricted to small distance quantities such as
the running coupling or scale dependent renomalization functions.

It would therefore be ideal to combine both approaches and perform quantum
computations at small values of the lattice spacing and large volume MCMC
simulations to provide the physical scale. To this end, in this paper we have

  * •

developed a resource efficient encoding for 2+1 dimensional QED, which allows
to simulate the model eventually on already now available or forthcoming near
term quantum computer with more and better qubits;

  * •

demonstrated that the introduction of suitable suppression terms can be used
to force the final state of the optimization of a generic ansatz to have a
small overlap with the unphysical Hilbert space;

  * •

obtained results for the ground state energy in broad range of coupling which
in turns allows us to compute the plaquette expectation value, or the static
force at small distances which can be related to the running coupling;

  * •

shown that we can obtain accurate enough results for the energy gap in the
large and intermediate coupling regime, which provides an important first step
to eventually to make the desired contact to MC simulations.

We also want to remark that the the setup of 2+1-dimensional QED developed
here will serve as the basis for extensions such as adding topological terms,
chemical potential or real time simulations, directions, we want to follow in
the future.

In this work we set the stage for an application of quantum computing (using
VQE in particular) to study lattice gauge theories in a non-perturbative
fashion. The present work addresses ground state properties as well as the
mass gap of 2+1-dimensional QED with the aim to reach small values of the
lattice spacing without running into problems with autocorrelations. However
the main advantage of this paradigm comes from the possibility of being
applied to study systems with numerical sign problems, which poses a challenge
to standard MCMC methods. In future works, given the here developed setup, we
will explore the addition of a chemical potential and topological $\theta$
term to the $(2+1)$-dimensional QED Hamiltonian, being thus able to go far
beyond traditional MCMC simulations – at least when quantum hardware will be
available to simulate large lattices.

###### Acknowledgements.

A.C. is supported in part by the Helmholtz Association - “Innopool Project
Variational Quantum Computer Simulations (VQCS)”. This work is supported with
funds from the Ministry of Science, Research and Culture of the State of
Brandenburg within the Centre for Quantum Technologies and Applications
(CQTA).

## References

  * [1] Y. Aoki et al., “FLAG Review 2021,” 11 2021.
  * [2] M. Troyer and U.-J. Wiese, “Computational complexity and fundamental limitations to fermionic quantum monte carlo simulations,” Phys. Rev. Lett., vol. 94, p. 170201, May 2005.
  * [3] M. A. Nielsen and I. Chuang, “Quantum computation and quantum information,” 2002\.
  * [4] O. Philipsen, “Lattice QCD at non-zero temperature and baryon density,” in Les Houches Summer School: Session 93: Modern perspectives in lattice QCD: Quantum field theory and high performance computing, pp. 273–330, 9 2010\.
  * [5] G. Aarts, “Introductory lectures on lattice QCD at nonzero baryon number,” J. Phys. Conf. Ser., vol. 706, no. 2, p. 022004, 2016.
  * [6] H.-T. Ding, F. Karsch, and S. Mukherjee, “Thermodynamics of strong-interaction matter from Lattice QCD,” Int. J. Mod. Phys. E, vol. 24, no. 10, p. 1530007, 2015.
  * [7] J. Preskill, “Quantum Computing in the NISQ era and beyond,” Quantum, vol. 2, p. 79, Aug. 2018.
  * [8] M. C. Bañuls et al., “Simulating Lattice Gauge Theories within Quantum Technologies,” Eur. Phys. J. D, vol. 74, no. 8, p. 165, 2020.
  * [9] T. Felser, P. Silvi, M. Collura, and S. Montangero, “Two-dimensional quantum-link lattice quantum electrodynamics at finite density,” Phys. Rev. X, vol. 10, p. 041040, Nov 2020.
  * [10] G. Mazzola, S. V. Mathis, G. Mazzola, and I. Tavernelli, “Gauge-invariant quantum circuits for $u$(1) and yang-mills lattice gauge theories,” Phys. Rev. Research, vol. 3, p. 043209, Dec 2021.
  * [11] E. Zohar, J. I. Cirac, and B. Reznik, “Simulating (2+ 1)-dimensional lattice qed with dynamical matter using ultracold atoms,” Physical review letters, vol. 110, no. 5, p. 055302, 2013.
  * [12] J. Ambjørn, A. J. Hey, and S. Otto, “String tensions for lattice gauge theories in 2+ 1 dimensions,” Nuclear Physics B, vol. 210, no. 3, pp. 347–368, 1982.
  * [13] M. Loan, M. Brunner, C. Sloggett, and C. Hamer, “Path integral monte carlo approach to the u (1) lattice gauge theory in 2+ 1 dimensions,” Physical Review D, vol. 68, no. 3, p. 034504, 2003.
  * [14] S. Booth, M. Göckeler, R. Horsley, A. Irving, B. Joo, S. Pickles, D. Pleiter, P. Rakow, G. Schierholz, Z. Sroczynski, and H. Stüben, “Determination of $\lambda_{MS}$ from quenched and nf=2 dynamical qcd,” Physics Letters B, vol. 519, no. 3, pp. 229–237, 2001.
  * [15] O. Raviv, Y. Shamir, and B. Svetitsky, “Nonperturbative beta function in three-dimensional electrodynamics,” Phys. Rev. D, vol. 90, p. 014512, Jul 2014.
  * [16] H. Nielsen and M. Ninomiya, “A no-go theorem for regularizing chiral fermions,” Physics Letters B, vol. 105, no. 2, pp. 219–223, 1981.
  * [17] L. Susskind, “Lattice fermions,” Phys. Rev. D, vol. 16, pp. 3031–3039, Nov 1977.
  * [18] H. J. Rothe, Lattice gauge theories: an introduction.  World Scientific Publishing Company, 2012.
  * [19] J. Kogut and L. Susskind, “Hamiltonian formulation of wilson’s lattice gauge theories,” Phys. Rev. D, vol. 11, pp. 395–408, Jan 1975.
  * [20] 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, vol. 5, p. 393, Feb. 2021.
  * [21] D. Paulson, L. Dellantonio, J. F. Haase, A. Celi, A. Kan, A. Jena, C. Kokail, R. van Bijnen, K. Jansen, P. Zoller, and C. A. Muschik, “Simulating 2d effects in lattice gauge theories on a quantum computer,” PRX Quantum, vol. 2, p. 030334, Aug 2021.
  * [22] K. G. Wilson, “Confinement of quarks,” Phys. Rev. D, vol. 10, pp. 2445–2459, Oct 1974.
  * [23] M. Bañuls, K. Cichy, J. Cirac, and K. Jansen, “The mass spectrum of the schwinger model with matrix product states,” Journal of High Energy Physics, vol. 2013, nov 2013.
  * [24] P. Jordan and E. P. Wigner, “Über das paulische äquivalenzverbot,” in The Collected Works of Eugene Paul Wigner, pp. 109–129, Springer, 1993\.
  * [25] O. Di Matteo, A. McCoy, P. Gysbers, T. Miyagi, R. M. Woloshyn, and P. Navrátil, “Improving hamiltonian encodings with the gray code,” Phys. Rev. A, vol. 103, p. 042405, Apr 2021.
  * [26] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, “A variational eigenvalue solver on a photonic quantum processor,” Nature Communications, vol. 5, jul 2014.
  * [27] G. Ortiz, J. E. Gubernatis, E. Knill, and R. Laflamme, “Quantum algorithms for fermionic simulations,” Physical Review A, vol. 64, no. 2, p. 022319, 2001\.
  * [28] J. Tilly et al., “The Variational Quantum Eigensolver: a review of methods and best practices,” 11 2021.
  * [29] T. Jones, S. Endo, S. McArdle, X. Yuan, and S. C. Benjamin, “Variational quantum algorithms for discovering hamiltonian spectra,” Phys. Rev. A, vol. 99, p. 062304, Jun 2019.
  * [30] O. Higgott, D. Wang, and S. Brierley, “Variational Quantum Computation of Excited States,” Quantum, vol. 3, p. 156, July 2019.
  * [31] V. Havlíček, A. D. Córcoles, K. Temme, A. W. Harrow, A. Kandala, J. M. Chow, and J. M. Gambetta, “Supervised learning with quantum-enhanced feature spaces,” Nature, vol. 567, pp. 209–212, mar 2019\.
  * [32] E. Knill, G. Ortiz, and R. D. Somma, “Optimal quantum measurements of expectation values of observables,” Phys. Rev. A, vol. 75, p. 012328, Jan 2007.
  * [33] K. M. Nakanishi, K. Fujii, and S. Todo, “Sequential minimal optimization for quantum-classical hybrid algorithms,” Phys. Rev. Research, vol. 2, p. 043158, Oct 2020.
  * [34] M. S. A. et al., “Qiskit: An open-source framework for quantum computing,” 2021\.
  * [35] R. Sommer, “A new way to set the energy scale in lattice gauge theories and its application to the static force and $\sigma$s in su (2) yang-mills theory,” Nuclear Physics B, vol. 411, no. 2-3, pp. 839–854, 1994.
  * [36] M. Loan, M. Brunner, C. Sloggett, and C. Hamer, “Path integral monte carlo approach to the u(1) lattice gauge theory in 2$+$1 dimensions,” Physical Review D, vol. 68, aug 2003.
  * [37] A. M. Polyakov, “Thermal Properties of Gauge Fields and Quark Liberation,” Phys. Lett. B, vol. 72, pp. 477–480, 1978.
  * [38] M. Göpfert and G. Mack, “Proof of Confinement of Static Quarks in Three-Dimensional U(1) Lattice Gauge Theory for All Values of the Coupling Constant,” Commun. Math. Phys., vol. 82, p. 545, 1981.
  * [39] M. Lüscher, P. Weisz, and U. Wolff, “A Numerical method to compute the running coupling in asymptotically free theories,” Nucl. Phys. B, vol. 359, pp. 221–243, 1991.
  * [40] M. Lüscher, R. Narayanan, P. Weisz, and U. Wolff, “The Schrödinger functional: A Renormalizable probe for nonAbelian gauge theories,” Nucl. Phys. B, vol. 384, pp. 168–228, 1992.
  * [41] M. Lüscher, “Advanced lattice QCD,” in Les Houches Summer School in Theoretical Physics, Session 68: Probing the Standard Model of Particle Interactions, pp. 229–280, 2 1998.

##  Appendix A Structured ansatz for incomplete Gray encoding

Here we discuss how a parameterized circuit should be adapted in order to
allow transitions only between encoded states.

For the cases where states are encoded using the computational basis, the
ansätze can be built by identifying the states to be excluded.

(a)  Figure 11: Parametric circuit for cycling between encoded states with the
Gray encoding at truncation $l=1$. (a)  Figure 12: Controlled parametric
circuit for cycling between encoded states with the Gray encoding at
truncation $l=1$.

For example, in the case of Gray encoding with $l=1$, shown in Table 1, one
can cycle between the three encoded states for each link
$\\{\ket{00},\ket{01},\ket{11}\\}$, using the parametric circuit shown in
Figure 11, where the two rotations parameterized by $\theta_{1}$ and
$\theta_{2}$ drive the transitions $\ket{00}\leftrightarrow\ket{01}$ and
$\ket{01}\leftrightarrow\ket{11}$ respectively, while the parametric circuit
shown in Figure 12 entangles the gauge register ($q_{1},q_{0}$) with a matter
component ($q_{c}$), by rotations which make for the transitions
$\ket{00}\leftrightarrow\ket{01}$ and $\ket{00}\leftrightarrow\ket{11}$ (the
rotation in $\theta_{2}$ is actually a $CR_{XX}(\theta_{2})$).

##  Appendix B Antipodal state as physical in the magnetic basis

As described in Sec. II.2.1, the discretization adopted in [20] to write the
Hamiltonian in the magnetic basis involves an uniform grid of $2L+1$ sites on
the group $U(1)$, where the identity is represented by
$\ket{0}_{\text{phys.}}$. From these states, only the $2l+1$ around the
identity are selected as physical ones. However, as already mentioned in Sec.
III.1, being the number of physical states not a power of two makes the
optimization problem harder, either at the Hamiltonian or at the ansatz level.
Including states non-symmetrically with respect to the identity (or origin in
the electric basis) would break the charge conjugation symmetry of the
Hamiltonian, therefore introducing another problem.

A straightforward solution, although working only in the magnetic case, would
be to change the discretization procedure: instead of using $2L+1$ sites
(i.e., the phases $\phi_{j}\equiv\frac{2\pi j}{2L+1}$), one can use a $2L+2$
grid, in such a way that the site antipodal to the identity is included in the
discretization. Then, at the truncation stage, besides the $2l+1$ states
around the identity we can include also the remaining
$2^{\lceil\log_{2}(2l+1)\rceil}-2l-1$ as states around the antipodal to the
identity. In this way, all the states used are encoded as physical states
belonging to the discretization grid and one can use a generic ansatz without
having to penalize unphysical states for the gauge variables.

Another advantage of this discretization into $2L+2$ states is that it is
possible to increase the truncation level from $L$ to $L^{\prime}=2L+1$ by
keeping the previous grid and inserting sites at its midpoints. Unlike the
$2L+1$ grid, which cannot be simply related to finer grids at larger $L$, this
would allow for a better control on the extrapolation to infinite
discretization.

##  Appendix C Penalty term method

We first discuss the introduction of a penalty term in the pure gauge system
with periodic system and then we mention some information about the case with
fermions.

###  C.1 Test in pure gauge theory

In the upper panel of Figure 13 it is depicted the trend of the energy
eigenvalues ($E_{0}$ in blue and $E_{1}$ in orange) for a certain range of
$\lambda$.

Figure 13: Suppression factor method for a given value of the coupling. Energy
eigenvalues from variational approach (optimization with 200 iterations and
average of 20 runs) and exact diagonalization (dotted lines in upper panel)
and amount of unphysical states in the VQD solution (lower panel).

When the suppression coefficient $\lambda$ is zero the first excited state is
almost degenerate to the ground level. However, one can see in the
corresponding value of unphysical states in the lower panel, that the result
does not represent a physical solution.

Figure 14: Suppression factor method for a given value of the coupling. With
800 iterations we increase the possibilities to obtain accurate results when
performing a variational method.

Figure 15: Trend of VQD results with fixed $\lambda=1000$ and different values
of iterations, (50, 100, 200 and 800).

As $\lambda$ increases, the percentage of unphysical states decreases, until
we can reach an accurate solution after the unphysical part of $E_{0}$ and
$E_{1}$ is suppressed, ($\lambda\sim O(10)$). However, it is important to not
have a large $\lambda$ when performing the variational method. This because
the algorithm cannot reach a low, acceptable value for the energies. In the
rightmost part of the plot the results deviate from the exact solution, even
if we can avoid unphysical states, the fidelity decreases888In some runs, it
may happen that the values found for the first excited is lower than the one
of the ground state. In this case one can exchange the two states and values..
If we consider more iterations for the optimization process, the possibilities
for the optimizer to converge increases even if $\lambda$ is large, as
depicted in Figure 14 and Figure 15.

Table 2: Energy eigenvalues $E_{k}$, in the coupling range $g\geq 1$ for electric basis with truncation $l=1$. Also the percentage of unphysical states, $u(\theta^{*}_{k})$ and number of layers $n_{\text{lay.}}$ in the quantum circuit are shown. $g$ | $n_{\text{lay.}}$ | $E_{0}$ | $E_{0}^{(\text{exact})}$ | $u(\theta^{*}_{0})$ | $E_{1}$ | $E_{1}^{(\text{exact})}$ | $u(\theta^{*}_{1})$
---|---|---|---|---|---|---|---
1.0 | 3 | -3.173(24) | -3.799 | 0.0 | -2.461(23) | -3.145 | 0.39
1.16993 | 4 | -2.272(24) | -2.915 | 0.195 | -1.410(18) | -2.075 | 0.293
1.36874 | 4 | -1.89(2) | -2.328 | 0.391 | -0.715(18) | -1.315 | 0.195
1.60133 | 5 | -1.59(2) | -1.911 | 0.098 | -0.073(18) | -0.561 | 0.0
1.87344 | 4 | -1.324(21) | -1.568 | 0.098 | 0.75(2) | 0.279 | 0.195
2.1918 | 5 | -1.085(23) | -1.265 | 0.0 | 1.539(27) | 1.251 | 0.098
2.56425 | 5 | -0.847(25) | -0.998 | 0.0 | 2.567(36) | 2.412 | 0.098
3.0 | 5 | -0.698(27) | -0.769 | 0.0 | 3.841(37) | 3.844 | 0.0

###  C.2 Penalty term for the fermionic system

For the fermionic case we applied a specific approach and considered only a
small range for the suppression parameter. In Table 2 we reported the best
results for the ground state energy and first excited state for the truncation
level $l=1$ and in the coupling range $g\in[1,3]$. The main criterion was to
set a threshold for the percentage of unphysical states. No explicit
dependence on the number of layers or small variations of $\lambda$ was found;
this might be due to the variability of the optimization process.