# First-Order Phase Transition of the Schwinger Model with a Quantum Computer

Takis Angelides  Institut für Physik, Humboldt-Universität zu Berlin,
Newtonstr. 15, 12489 Berlin, Germany  CQTA, Deutsches Elektronen-Synchrotron
DESY, Platanenallee 6, 15738 Zeuthen, Germany  Pranay Naredi  Computation-
Based Science and Technology Research Center, The Cyprus Institute, 20 Kavafi
Street, 2121 Nicosia, Cyprus  Arianna Crippa  Institut für Physik, Humboldt-
Universität zu Berlin, Newtonstr. 15, 12489 Berlin, Germany  CQTA, Deutsches
Elektronen-Synchrotron DESY, Platanenallee 6, 15738 Zeuthen, Germany  Karl
Jansen  CQTA, Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, 15738
Zeuthen, Germany  Computation-Based Science and Technology Research Center,
The Cyprus Institute, 20 Kavafi Street, 2121 Nicosia, Cyprus  Stefan Kühn
CQTA, Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, 15738 Zeuthen,
Germany  Ivano Tavernelli  IBM Research Europe - Zurich, Säumerstrasse 4, 8803
Rüschlikon, Switzerland  Derek S. Wang  IBM Quantum, IBM T.J. Watson Research
Center, Yorktown Heights, NY 10598, USA

###### Abstract

We explore the first-order phase transition in the lattice Schwinger model in
the presence of a topological $\theta$-term by means of the variational
quantum eigensolver (VQE). Using two different fermion discretizations, Wilson
and staggered fermions, we develop parametric ansatz circuits suitable for
both discretizations, and compare their performance by simulating classically
an ideal VQE optimization in the absence of noise. The states obtained by the
classical simulation are then prepared on the IBM’s superconducting quantum
hardware. Applying state-of-the art error-mitigation methods, we show that the
electric field density and particle number, observables which reveal the phase
structure of the model, can be reliably obtained from the quantum hardware. To
investigate the minimum system sizes required for a continuum extrapolation,
we study the continuum limit using matrix product states, and compare our
results to continuum mass perturbation theory. We demonstrate that taking the
additive mass renormalization into account is vital for enhancing the
precision that can be obtained with smaller system sizes. Furthermore, for the
observables we investigate we observe universality, and both fermion
discretizations produce the same continuum limit.

##  I Introduction

Exploring the phase diagram of a theory is central for the understanding of
the fundamental laws of physics in many disciplines ranging from condensed
matter to particle physics. For instance, in the case of ferromagnetic
materials, critical points and phase boundaries help to understand and design
magnetic materials with various applications [1, 2]. Within superconductivity,
phase diagrams allow for identifying the critical temperatures and magnetic
fields required for materials to exhibit zero electrical resistance [3, 4]. In
particle physics, the phase diagram of Quantum Chromodynamics (QCD) elucidates
the behavior of matter under extreme conditions, relevant for the early
universe and neutron stars [5]. In high-energy physics, topological terms play
an important role in the study of interesting phenomena such as the breaking
of charge conjugation-parity symmetry in QCD [6], and the occurrence of out-
of-equilibrium dynamical effects involving axion fields [7].

In order to explore the phase structure of a given theory, one often has to
resort to numerical methods. A standard tool is Monte Carlo (MC) simulation,
which has been successful for probing phase diagrams [8, 9, 10, 11], including
that of lattice QCD [12]. However, the sign problem poses a barrier in certain
parameter regimes [13], leaving relevant questions unanswered. Specifically
for lattice QCD, large baryon chemical potentials or the presence of a
topological $\theta$-term would trigger this problem [14]. Thus, variations of
the conventional MC approach have been proposed to tackle this obstacle [12],
but so far only with limited success. Hence, there is great interest in
alternative methods that bypass the issue, such as tensor networks and quantum
computing [15, 16, 17]. In particular, quantum computing offers a promising
alternative route [16, 17, 18], with already a number of successful
demonstrations [19, 20, 21, 22, 23, 24].

Here, we focus on exploring the possibility of studying the phase structure of
the Schwinger model [25] in the presence of a topological $\theta$ term with
near-term quantum devices. Despite its simplicity, the Schwinger model in the
presence of a topological $\theta$-term provides an example where conventional
Monte Carlo methods would suffer from the sign problem. Moreover, it exhibits
a rich phase structure with a first-order quantum phase transition for large
enough masses at $\theta=\pi$ [26, 27, 28, 29, 30, 31, 32, 33, 34]. Hence,
Schwinger model serves as a benchmark system for developing and testing new
methods. Regarding a quantum computing approach to the model in the presence
of a topological $\theta$-term, there are several open questions. First, the
theory needs to be discretized on a lattice with a finite extent, and there
are various different ways for discretizing fermions on such a lattice. It is
therefore a priori not clear which one will show the best performance for a
given set of resources. Second, it is important to identify how large the
system sizes should be in order to obtain a reliable continuum extrapolation
in order to assess the applicability of near-term quantum devices.

Here we address these questions by studying two fermion discretizations,
staggered and Wilson fermions, in order to explore any potential advantages of
one type of fermion formulation over the other. We use the variational quantum
eigensolver (VQE) [35] and derive a protocol for mapping this problem on a
quantum circuit. The VQE, introduced as an alternative to quantum phase
estimation [36], aligns with the capabilities of current and near-term quantum
devices. We test our VQE using noiseless classical simulations between 6 to 12
qubits, to identify the best possible setup regarding the ansatz and gates
that would capture the relevant ground states most efficiently. After the
optimal ansatz-gate combination and variational parameters are found, the
ground states across the phase transition are prepared on IBM’s quantum
devices. We deomstrate that using state-of-the-art error mitigation
techniques—zero noise extrapolation [37], readout error mitigation [38], Pauli
twirling [39] and dynamical decoupling [40]—allows for obtaining precise
results from the quantum measurements. To understand the minimum system sizes
required to extrapolate faithfully to the continuum limit with a quantum
computer, we use matrix product states (MPS). We numerically simulate
intermediate system sizes and perform the continuum extrapolation, which we
compare to analytical results. Our study also shows universality for the
considered observables, as both discretizations lead to the same continuum
values.

The rest of the paper is structured as follows. In Sec. II, we briefly
introduce the Schwinger model and review its phase structure in the presence
of a topological $\theta$-term. Moreover, we discuss two different ways of
discretizing it on a lattice using Wilson and staggered fermions. We proceed
with presenting our ansatz for the VQE as well as the MPS techniques we use to
estimate the resources for taking a reliable continuum limit in Sec. III. Our
numerical results demonstrating the performance of the ansatz in various
parameter regimes are presented in Sec. IV, before concluding in Sec. V.

##  II The Schwinger model

The Schwinger model describes quantum electrodynamics in (1+1)-dimensions
coupled to a single, massive Dirac fermion [25]. Here we briefly introduce the
Hamiltonian formulation and review its phase diagram in the presence of a
topological $\theta$-term. We then discuss two different discretizations for
the fermionic matter fields of the model, namely Wilson and staggered
fermions.

###  II.1 Hamiltonian formulation in the continuum

The continuum Hamiltonian density of the Schwinger model in the presence of a
topological $\theta$-term is given by

$\mathcal{H}=-i\overline{\psi}\gamma^{1}\left(\partial_{1}-igA_{1}\right)\psi+m\overline{\psi}\psi+\frac{1}{2}\left(\dot{A}_{1}+\frac{g\theta}{2\pi}\right)^{2},$
(1)

where $\psi(x)$ is a two-component Dirac spinor describing the fermionic
matter. The spinor components, $\psi_{\alpha}$, $\alpha=1,2$, fulfill the
standard fermionic anticommutation relations
$\\{\psi_{\alpha}^{\dagger}(x),\psi_{\beta}(y)\\}=\delta(x-y)\delta_{\alpha\beta}$.
The gauge field $A_{\mu}$, $\mu=0,1$, mediates the interaction between the
matter fields. Here we have chosen the temporal gauge, $A_{0}=0$, hence only
the spatial component $A_{1}$ appears in the Hamiltonian. The parameters $m$
and $g$ are the bare fermion mass and the coupling between fermions and the
gauge fields. The matrices $\gamma^{\mu}$ are two dimensional matrices obeying
the Clifford algebra $\\{\gamma^{\mu},\gamma^{\nu}\\}=2\eta^{\mu\nu}$, with
$\eta=\text{diag}(1,-1)$, and $\overline{\psi}$ corresponds to
$\psi^{\dagger}\gamma^{0}$. The physically relevant gauge invariant states of
the Hamiltonian have to fulfill Gauss’s law

$\displaystyle-\partial_{1}\dot{A}^{1}=g\overline{\psi}\gamma^{0}\psi,$ (2)

where $-\dot{A}^{1}$ is the electric field and
$g\overline{\psi}\gamma^{0}\psi$ represents the charge density.

The topological term, $g\theta/2\pi$, appearing in the Hamiltonian corresponds
to a constant background electric field whose effect has been assessed both
theoretically and numerically. Coleman argued that the physics of the model is
periodic in $\theta$ with a period of $2\pi$, and that above a certain
critical mass, $m_{c}/g$, the model undergoes a first order quantum phase
transition at $\theta=\pi$ [30]. This picture was later on confirmed in
numerical simulations, where it was found that the critical line ends in a
second-order quantum phase transition at $m_{c}/g\sim 0.33$ [29, 34, 41].
Figure 1 provides a sketch of the phase diagram, highlighting the first-order
phase transition line, which culminates with a second-order phase transition
at $m_{c}/g$.

Figure 1: Illustration of the phase diagram of the Schwinger model in the
presence of a topological term in the $m/g-\theta$ plane. Since the physics is
periodic in $\theta$ with period $2\pi$, only the first period is shown. The
critical line (shown in black) indicates the first-order phase transitions
occurring at $\theta=\pi$ for masses larger than the critical one
$m_{c}/g\approx 0.33$, which ends in a second-order phase transition (green
dot) exactly at $m_{c}/g$. Below the critical mass no transitions occur.

The physics of the model can also be understood qualitatively in an intuitive
picture. For large values of the mass in units of the coupling, $m/g\gg 1$, it
is generally unfavorable to generate charged particles. In the regime
$l_{0}=\theta/2\pi<1/2$, the ground state thus remains devoid of charged
particles and the electric field matches the background electric field due to
Gauss’s law. For $l_{0}>1/2$, it becomes energetically favorable to produce a
negative unit charge on the left of the one-dimensional space and a particle
of positive unit charge on the right. As a result of Gauss’s law, these are
connected by a flux string, lowering the electric field by one unit, which
takes the electric field in the bulk from $l_{0}$ to $l_{0}-1$ [34, 30]. Given
that the electric field energy is proportional to the square of the electric
field and the mass contribution of the particle-antiparticle pair is
negligible at infinite volume, these two states are degenerate directly at the
point $l_{0}=1/2$. Hence, a first-order quantum phase transition occurs at
$l_{0}=1/2$ (or equivalently $\theta=\pi$) for large values of $m/g$,
manifesting itself in an abrupt change of the electric field density and
particle number. This transition spontaneously breaks the invariance under a
combined charge conjugation and parity symmetry transformation (CP
transformation) [41].

For small values of $m/g$ below the critical mass, it is energetically
favorable to screen the electric field by producing more of the negative-
positive charge pairs, and at $m/g=0$ we have complete screening of the
background electric field [42]. Hence, in this regime the electric field
density and the particle number do not display an abrupt change as a function
of $\theta$.

###  II.2 Lattice discretizations

In order to address the Hamiltonian of the theory numerically, we choose to
work with a discrete lattice formulation. To ensure that one faithfully
recovers the physics of the continuum model in the limit of vanishing lattice
spacing, special care has to be taken when discretizing the fermionic degrees
of freedom. In particular, a naive discretization of the fermionic fields
leads to an incorrect continuum limit due to the so-called doubling problem
[43]. There are several ways to avoid this problems, and here we choose to
work with two types of fermion discretizations, Wilson and staggered fermions,
which are commonly used in the literature. We first turn to Wilson fermions
and present the final Hamiltonian as well as the relevant observables in this
formulation before discussing the same operators for the staggered
discretization. Moreover, we show how to translate them to qubit degrees of
freedom.

####  II.2.1 Wilson fermions

The basic idea of Wilson fermions is to add a lattice version of the second
derivative of the fermion field to the native discretization (for details, see
Appendix A). While this term vanishes in the continuum limit, it can be shown
that it breaks the chiral symmetry explicitly, and gives the doubler solutions
originating from the discretization of space a mass proportional to the
inverse of the lattice spacing. This ensures that they decouple from the
theory when approaching the continuum limit [44, 45]. The Wilson
discretization of the of the Schwinger Hamiltonian on a lattice with $N$ sites
and spacing $a$ reads [26, 46]

$\displaystyle\begin{aligned}
H_{W}&=\sum_{n=0}^{N-2}\left(\overline{\phi}_{n}\left(\frac{r+i\gamma^{1}}{2a}\right)U_{n}\phi_{n+1}+\text{h.c.}\right)\\\
&+\sum_{n=0}^{N-1}\left(m_{\text{lat}}+\frac{r}{a}\right)\overline{\phi}_{n}\phi_{n}+\sum_{n=0}^{N-2}\frac{ag^{2}}{2}\left(L_{n}+l_{0}\right)^{2}.\end{aligned}$
(3)

The field $\phi_{n}$ is a dimensionless two-component Dirac spinor at site
$n$, whose components $\phi_{n,\alpha}$, $\alpha=1,2$, again fulfill the
standard anticommutation relations for fermionc fields
$\\{\phi_{n,\alpha}^{\dagger},\phi_{n^{\prime},\beta}\\}=\delta_{nn^{\prime}}\delta_{\alpha\beta}$.
The operators $U_{n}$ and $L_{n}$ act on the links in between the matter sites
$n$ and $n+1$, where $L_{n}$ represents the quantized dimensionless electric
field. They fulfill the commutation relation
$[U_{n},L_{n^{\prime}}]=\delta_{nn^{\prime}}U_{n}$, hence, in the eigenbasis
of the electric field operator, $U_{n}$ acts as a lowering operator,
decreasing the electric flux value by one unit. The parameters
$m_{\text{lat}}$ and $g$ are the lattice mass and coupling, and $r$
corresponds the Wilson parameter, which can be chosen arbitrarily between
$(0,1]$ to ensure the correct continuum limit [46]. The background electric
field, $l_{0}=\theta/2\pi$, is the lattice version of the topological
$\theta$-term. On the lattice, Gauss’s law translates to

$\displaystyle L_{n}-L_{n-1}=Q_{n}$ (4)

where $Q_{n}=\phi_{n}^{\dagger}\phi_{n}-1$ is the charge operator. Given the
global $U(1)$ symmetry of the model, the total charge $\sum_{n}Q_{n}$ is
conserved. The lattice discretization with Wilson formulation generates an
additive mass renormalization [26]. The renormalized mass can be expressed as

$m_{r}/g=m_{\text{lat}}/g+\text{MS}(V,ag,l_{0}),$ (5)

where MS is the additive mass shift which, in general, depends on the
dimensionless lattice volume, $V=N/\sqrt{x}$, the lattice spacing in units of
the coupling, $ag$, and the background field, $l_{0}$. In order to follow a
line of constant $m_{\text{r}}/g$, one needs to determine MS. For our
simulations later on, we use the method from Ref. [26] to obtain MS. For the
Wilson discretization, we can define the particle number operator as

$P_{W}=N+\sum_{n=0}^{N-1}\overline{\phi}_{n}\phi_{n},$ (6)

which essentially is the lattice version of the chiral condensate,
$\overline{\psi}\psi$, up to the constant $N$. The constant shift ensures that
$P_{W}$ is a positive semidefinite operator, where the zero eigenvalue is
obtained for states minimizing $\sum_{n=0}^{N-1}\overline{\phi}_{n}\phi_{n}$.
These correspond to the ground state of the Hamiltonian in the limit of
dominating mass term, for which we do not expect any particles from the
picture discussed for the continuum model. As soon as particles are present in
the ground state, we expect a nonzero expectation value for $P_{W}$. Thus,
together with the electric field, this observable should allow us to detect
the first-order phase transition in the lattice model.

For open boundary conditions, Eq. (4) can be solved iteratively after fixing
the value of the electric field on the left boundary, $\varepsilon_{0}$.
Choosing $\varepsilon_{0}=0$, we find for the electric field111Choosing a
nonzero $\varepsilon_{0}$ value for the electric field can be interpreted as a
shift in the value of $l_{0}$ to $l_{0}^{\prime}=l_{0}+\varepsilon_{0}$.
Hence, we can set $\varepsilon_{0}=0$ without loss of generality.

$\displaystyle L_{n}=\sum_{k=0}^{n}Q_{k}.$ (7)

This implies that the electric field on each link can be determined by the
fermionic charges, and is not a dynamical variable. Inserting this into the
Hamiltonian in Eq. (3) and applying a residual transformation to the fermionic
fields [47], we can integrate out the gauge fields and obtain a formulation
purely in fermionic degrees of freedom [26]

$\displaystyle\begin{aligned}
H_{W}&=\sum_{n=0}^{N-2}\left(\overline{\phi}_{n}\left(\frac{r+i\gamma^{1}}{2a}\right)\phi_{n+1}+\text{h.c.}\right)\\\
&+\sum_{n=0}^{N-1}\left(m_{\text{lat}}+\frac{r}{a}\right)\overline{\phi}_{n}\phi_{n}+\sum_{n=0}^{N-2}\frac{ag^{2}}{2}\left(l_{0}+\sum_{k=0}^{n}Q_{k}\right)^{2}.\end{aligned}$
(8)

Note that the above Hamiltonian is directly confined to the gauge invariant
subspace, as Gauss’s law is fulfilled by construction. For the rest of this
paper, we will fix $r=1$ and choose to work with the convention
$\gamma^{0}=X$, $\gamma^{1}=iZ$, where $X$, $Z$ are the usual Pauli matrices.

In order to be able to measure the Hamiltonian on a quantum device, we have to
map the fermionic degrees of freedom to qubits. To this end, we use the
Jordan-Wigner transformation [48] to translate the fermion fields to spins.
Note that in the Wilson formulation, we have a two-component spinor on each
site. Hence, the Jordan-Wigner transformation for a system with $N$ sites
results in $2N$ spin degrees of freedom. As shown in detail in Appendix A,
applying a convenient ordering for the fermions, allows for obtaining a purely
real Hamiltonian, and we find for the dimensionless spin formulation of Eq.
(8)

$\displaystyle\begin{aligned}
W_{W}&=x\sum_{n=0}^{N-2}\left(X_{2n+2}X_{2n+3}+Y_{2n+2}Y_{2n+3}\right)\\\
&+\left(\frac{m_{\text{lat}}}{g}\sqrt{x}+x\right)\sum_{n=0}^{N-1}\left(X_{2n+1}X_{2n+2}+Y_{2n+1}Y_{2n+2}\right)\\\
&+\sum_{n=0}^{N-2}\left(l_{0}+\sum_{k=0}^{n}Q_{k}\right)^{2}.\end{aligned}$
(9)

Here we have defined the inverse lattice spacing squared in units of the
coupling, $x=1/(ag)^{2}$. The charge operator in spin formulation is given by
$Q_{n}=(Z_{2n}+Z_{2n+1})/2$, and $X_{n}$, $Y_{n}$, and $Z_{n}$ are the
standard Pauli operators acting on spin $n$.

The relevant observables for detecting the first-order phase transition, can
be readily translated to the spin formulation. Since we expect the electric
field density to be homogeneous in the bulk for large enough system sizes, we
only consider a single link operator in the center of the system to minimize
boundary effects. The total electric field including the background field on
this link can then be reconstructed from the charges as

$\displaystyle L_{W}$
$\displaystyle=l_{0}+L_{\lceil\frac{N}{2}\rceil-1}=l_{0}+\sum_{k=0}^{\lceil
N/2\rceil-1}Q_{k}$ (10) $\displaystyle=l_{0}+\frac{1}{2}\sum_{k=0}^{\lceil
N/2\rceil-1}(Z_{2k}+Z_{2k+1}).$

Using the same ordering for the fermions as for the Hamiltonian, the particle
number operator becomes

$P_{W}=N+\frac{1}{2}\sum_{n=0}^{N-1}\left(X_{2n+1}X_{2n+2}+Y_{2n+1}Y_{2n+2}\right).$
(11)

This operator commutes with the Hamiltonian only in the weak-coupling limit,
$m/g\gg 1$, for which the kinetic term can be neglected. Hence, only in this
limit we expect the integer eigenvalues of $P_{W}$ to be good quantum numbers
[49].

####  II.2.2 Staggered fermions

Another approach to discretize the model are the Kogut-Susskind staggered
fermions [50]. The basic idea is to distribute the two spinor components to
different lattice sites, thereby “thinning out” the fermionic degrees of
freedom and avoiding the doubling problem. The Kogut-Susskind formulation for
the Schwinger Hamiltonian for a lattice with $N$ sites, spacing $a$ and open
boundary conditions can written as [50]

$\displaystyle\begin{aligned}
H_{S}&=-\frac{i}{2a}\sum_{n=0}^{N-2}\left(\phi^{\dagger}_{n}U_{n}\phi_{n+1}-\text{h.c}.\right)\\\
&+m_{\text{lat}}\sum_{n=0}^{N-1}(-1)^{n}\phi^{\dagger}_{n}\phi_{n}+\frac{ag^{2}}{2}\sum_{n=0}^{N-2}L_{n}^{2}\,,\end{aligned}$
(12)

where $N$ is supposed to be even. In the expression above, the operators
$L_{n}$ and $U_{n}$ are the same as for the Wilson case and act on the links
between sites $n$ and $n+1$. Different to the Wilson discretization,
$\phi_{n}$ now represents a single-component fermionic field whose
anticommutation relations are given by
$\\{\phi_{n}^{\dagger},\phi_{n^{\prime}}\\}=\delta_{nn^{\prime}}$. In the
limit of vanishing lattice spacing, the fermionic fields residing on even
(odd) sites will correspond to the upper (lower) component of the Dirac spinor
$\psi(x)$. Thus, the staggered formulation effectively doubles the lattice
spacing. Again, $m_{\text{lat}}$ and $g$ refer to the lattice mass and the
coupling. The Gauss’s law for the staggered case has the same form as in Eq.
(4), but $Q_{n}$ is now given by the staggered charge operator
$Q_{n}=\phi^{\dagger}_{n}\phi_{n}-\left(1-(-1)^{n}\right)/2$. Also for the
staggered formulation, the total charge $\sum_{n}Q_{n}$ is a conserved
quantity and commutes with the Hamiltonian.

Similar to the Wilson case, it was recently shown that staggered fermions
experience an additive mass renormalization, where the renormalized mass can
again be written as in Eq. (5) [51]. For the case of periodic boundary
conditions, Ref. [51] was able to compute the mass shift analytically finding

$\displaystyle\text{MS}_{t}=\frac{1}{8\sqrt{x}}.$ (13)

For open boundary conditions, there is no analytical prediction, and we can
measure the mass shift using the approach from Ref. [26]. Moreover, Ref. [26]
also demonstrated that the results for MS obtained for the staggered
discretization using open boundary conditions agree with the theoretical
prediction $\text{MS}_{t}$ for large lattice volumes. Thus, for large enough
volumes Eq. (13) should still provide an approximation even for open boundary
conditions. In our numerical computations later on, we will consider both
options.

Analogously to the Wilson case, the particle number operator for the staggered
fermions resembles the mass term of the Hamiltonian and is given by

$P_{S}=\frac{N}{2}+\sum_{n=0}^{N-1}(-1)^{n}\phi_{n}^{\dagger}\phi_{n}$ (14)

This expression is again equivalent the lattice version of the chiral
condensate up to a constant shift, and has the same properties as for the
Wilson case. Note that $P_{S}$ now contains an alternating sign factor, which
is the result of the staggered discretization separating the spinor components
to different lattice sites. Together with the electric field, this observable
is again suitable to detect the first-order phase transition in the model.

As in the case of Wilson fermions, the gauge fields can be integrated using
Gauss’s law resulting in Eq. (7), but now with the staggered charge.
Performing the same steps as in the Wilson case, we can obtain a Hamiltonian
directly restricted to the gauge invariant subspace where all notion of the
gauge field is gone [47]

$\displaystyle\begin{aligned}
H_{S}&=-\frac{i}{2a}\sum_{n=0}^{N-2}\left(\phi^{\dagger}_{n}\phi_{n+1}-\text{h.c}.\right)\\\
&+m_{\text{lat}}\sum_{n=0}^{N-1}(-1)^{n}\phi^{\dagger}_{n}\phi_{n}+\frac{ag^{2}}{2}\sum_{n=0}^{N-2}\left(l_{0}+\sum_{k=0}^{n}Q_{k}\right)^{2}.\end{aligned}$
(15)

Finally, this expression can again be translated to qubits by using the
Jordan-Wigner transformation [48] (for details, see Appendix A). Note that in
the staggered case, each site is only populated by a single-component
fermionic field. Thus, in contrast to the Wilson case, for $N$ physical sites
we obtain a system with $N$ qubits and the dimensionless Hamiltonian reads
[47]

$\displaystyle\begin{aligned}
W_{S}&=\frac{x}{2}\sum_{n=0}^{N-2}\left(X_{n}X_{n+1}+Y_{n}Y_{n+1}\right)\\\
&+\frac{m_{\text{lat}}}{g}\sqrt{x}\sum_{n=0}^{N-1}(-1)^{n}Z_{n}+\sum_{n=0}^{N-2}\left(l_{0}+\sum_{k=0}^{n}Q_{k}\right)^{2}\end{aligned}$
(16)

where $x=1/(ag)^{2}$ is again the inverse lattice spacing in units of the
coupling squared, and $X_{n}$, $Y_{n}$, and $Z_{n}$ are the Pauli matrices
acting on spin $n$ as in the Wilson case. The staggered charge operator is
given by $Q_{n}=(Z_{n}+(-1)^{n})/2$, from which we can infer that on even
(odd) sites we have a positron (electron) present when the spin is up (down).

The observables for detecting the first-order phase transition can again be
straightforwardly expressed in spin language. As opposed to Wilson fermions,
in the staggered discretization, the charge configuration on the lattice
exhibits a staggering effect, which leads to a non-uniform electric flux on
the links, in particular for smaller lattices. Therefore, to suppress the
boundary effects as well as the staggering effect while computing the electric
field density, two adjacent links in the center of the system are averaged
rather than just a single one for the case of Wilson fermions. The
corresponding electric field operator including the background field is given
by

$\displaystyle\begin{split}L_{S}&=l_{0}+\frac{1}{2}\left(L_{N/2-2}+L_{N/2-1}\right)\\\
&=l_{0}+\frac{1}{2}\left(\sum_{k=0}^{N/2-2}Q_{k}+\sum_{k=0}^{N/2-1}Q_{k}\right)\\\
&=l_{0}+\frac{1}{4}+\frac{1}{2}\sum_{k=0}^{N/2-2}Z_{k}+\frac{1}{4}Z_{N/2-1}.\end{split}$
(17)

For the particle number we find

$\displaystyle P_{S}=\frac{N}{2}+\frac{1}{2}\sum_{n=0}^{N-1}(-1)^{n}Z_{n}.$
(18)

When there are no particles present, i.e. the state of the system is given by
spin down on even sites and spin up on odd sites, the sum in the equation
above contributes a $-1$ for every site. In contrast, for even sites with spin
up and on odd sites with spin down, the sum contributes a $+1$. Hence, the
first term in Eq. (18) is added to render $P_{S}$ positive semidefinite.

In summary, both staggered and Wilson fermions provide a viable discretization
of the continuum Schwinger Hamiltonian. For a lattice system with $N$ physical
sites the Wilson discretization results in a Hamiltonian on $2N$ qubits,
whereas the staggered approach only requires $N$ qubits. While both
formulations reproduce the correct continuum limit, as we demonstrate
explicitly in Sec. IV, it is a priori not clear which of the two
discretizations converges faster. Although the staggered approach requires
less qubits than the Wilson one, this does not imply it will produce better
results given a fixed amount of resources. The latter question is particularly
relevant for quantum computing, as current and near-term devices only offer a
limited number of qubits that still suffer from a considerable level of noise.
Here, we aim at testing both approaches in a realistic scenario to benchmark
their performance. For the rest of the paper we will focus on the sector of
vanishing total charge, $\sum_{n}Q_{n}=0$, for both approaches.

##  III Methods

In order to assess both fermion discretizations, we study their performance
with a VQE as well as their convergence towards the continuum limit. Here, we
introduce the VQE setup we consider, including a description of the parametric
ansatz circuits and the optimization procedure for the parameters we utilize.
Moreover, we discuss the error mitigation techniques used for the inference
runs on quantum hardware. Finally, we briefly describe the MPS techniques we
use to explore the behavior of both discretizations towards the continuum
limit.

###  III.1 Parametric ansatz circuits for VQE

In order to test the performance of the different discretizations for VQE, we
focus on two different types of parametric ansätze and consider two different
types of gates, as shown in Fig. 2. We refer to the two ansatz architectures
as “brick” (c.f. Fig. 2(c)) and “ladder” (c.f. Fig. 2(d)). The two types of
parametric gates we consider are the SO(4) gates and the $R_{XX+YY}$ gates,
whose decomposition in standard controlled-NOT (CNOT) and Pauli rotation gates
is shown in Figs. 2(a) and 2(b). Here, we have chosen SO(4) instead of SU(4)
gates, because the Hamiltonians we study are real, and hence their ground
states are real. Thus, we can restrict our ansätze to the real subspace of the
Hilbert space. While the SO(4) gates are in principle more expressive than the
$R_{XX+YY}$ ones, they do not conserve the total charge. As a result, if we
use SO(4) gates, we need to manually enforce vanishing total charge, which we
do by adding a penalty term $\lambda\left(\sum_{n=0}^{N-1}Q_{n}\right)^{2}$ to
the Hamiltonians in Eq. (9) and Eq. (16). The Lagrange multiplier $\lambda$
has to be chosen sufficiently large that one obtains a ground state with
vanishing total charge. In contrast, the $R_{XX+YY}$ gate preserves the total
charge, but is generally less expressive.

Figure 2: Decomposition of a generic SO(4) gate depending on the six
parameters $\theta_{1},\dots,\theta_{6}$ (a) the
$R_{XX+YY}(\theta)=\exp(-i\theta(XX+YY)/2)$ (b), into CNOT and Pauli rotation
gates. Boxes acting on a single qubit correspond to Pauli rotation gates,
$R_{P}(\alpha)=\exp{-i\alpha P/2}$ with $P\in\\{X,Y,Z\\}$. Single-qubit gates
where the argument is omitted refer to rotations around an angle $\pi/2$,
$R_{P}(\pi/2)$. The light blue boxes represent the parameterized gates which
are $R(\alpha,\beta,\gamma)=R_{X}(\gamma)R_{Z}(\beta)R_{X}(\alpha)$ in (a) and
$R_{Y}(\theta)$ in (b). Panel (c) and (d) illustrate one layer of the brick
and ladder ansatz, respectively, both following a non-parametric part for
preparing the initial state $\ket{\psi_{in}}$ (yellow box). The first layer in
the brick ansatz has a CNOT-depth of $4$ whereas in ladder it is $2n-2$, where
$n$ is the number of qubits, and in both cases it increases by $4$ with each
layer.

In all our simulations, we choose the initial parameters for the parametric
part of this ansatz randomly in the interval $[0,0.001)$, such that the
parametric part of the circuit is close to the identity and we start the VQE
with a state close to $\ket{\psi_{in}}$.

For Wilson fermions, three options for the ansatz’s non-parametric portion
$\ket{\psi_{in}}$ are tested, and the one producing the best fidelity in the
VQE will be adopted. As explained in Sec. II.1, for $m/g\gg 1$ and
$\theta<\pi$, we expect the continuum Hamiltonian to be dominated by the mass
term, and the ground state corresponds to a state with no particles present.
The lattice analog of this state for the Wilson case is given by the product
state $\left(({\ket{01}-\ket{10}})/\sqrt{2}\right)^{\otimes N}$. As one can
easily verify, this state contains no charges and should provide a good
initial state for large masses before the onset of the first-order phase
transition. For $\theta>\pi$ the continuum model predicts that a particle
antiparticle pair is produced and a negative (positive) charge will form on
the left (right) boundary of the system. On the lattice, this would correspond
to the state
$\ket{11}\left(({\ket{01}-\ket{10}})/\sqrt{2}\right)^{\otimes(N-2)}\ket{00}$
in spin formulation, which approximates the state expected after the phase
transition. Finally, we also try a generic state with vanishing total charge
given by $\ket{10}^{\otimes N}$. In all cases, a simple linear mapping of the
logical qubits $\ket{\sigma_{l}}$ in the model Hamiltonian to the physical
ones of the hardware $\ket{q_{l}}$ is used.

For the staggered fermions, the state $\ket{10}^{\otimes N/2}$ corresponds to
the ground state at large values of $m/g\gg 1$ in the sector of vanishing
total charge and, thus, is a viable initial state $\ket{\psi_{in}}$ for large
masses and $\theta<\pi$. For $\theta>\pi$ we expect similar to the Wilson case
a pair of charges to form with the negative (positive) charge located on the
left (right) boundary of the system. Due to the staggered formulation,
negative (positive) charges can only reside on odd (even) sites, thus the pair
will be located at sites $1$ and $N-2$, resulting in the state
$\ket{11}\ket{10}^{\otimes(N-2)/2}\ket{00}$. This state can be easily created
from $\ket{10}^{\otimes N/2}$ by acting with an SO(4) or an $R_{XX+YY}$ gate
between qubits $1$ and $N-2$. For the ansatz architectures shown in Fig. 2,
there is no two-qubit gate directly acting on qubits $1$ and $N-2$. In order
to provide such a connection between these qubits, we deviate from the usual
linear mapping of the logical qubits in the model Hamiltonian to the physical
ones of the hardware. In particular, we choose the mapping
$\ket{q_{0}}\otimes\ket{q_{1}}\otimes\dots\otimes\ket{q_{N-1}}\equiv\ket{\sigma_{(N-2)/2+1}}\otimes\dots\otimes\ket{\sigma_{N-3}}\otimes\ket{\sigma_{N-1}}\otimes\ket{\sigma_{N-2}}\otimes\ket{\sigma_{1}}\otimes\ket{\sigma_{0}}\otimes\ket{\sigma_{2}}\otimes\dots\otimes\ket{\sigma_{(N-2)/2}}$,
where $\ket{q_{k}}$, $q_{k}\in\\{0,1\\}$ are the basis states of the qubits
$k$, and $\ket{\sigma_{l}}$, $\sigma_{l}\in\\{0,1\\}$ represent the basis
states of the logical qubit $l$ in the Hamiltonian. Using the ansätze shown in
Fig. 2 on the physical qubits, this provides a direct two-qubit gate between
the logical qubits $1$ and $N-2$. Thus, for the staggered case with large
masses, we use the initial state $\ket{10}^{\otimes N/2}$ for the logical
qubits and then proceed with the mapping discussed above. In the opposite
regime, $m/g\ll 1$, in which we expect the ground state to (partially) screen
the electric field, we do not foresee that a deviation from the linear mapping
is required. Thus in this parameter range, we choose $\ket{10}^{\otimes N/2}$
as an initial state for the logical qubits and simply use the linear mapping
as in the Wilson case.

In order to benchmark the ansätze for the different discretizations, we first
simulate the VQE classically. We employ the L-BFGS-B optimizer [52] with two
distinct heuristic warm-start stages. In the first phase, we group all
variational parameters in a layer and collectively optimize them using a VQE
with 2000 iterations. This stage is based on the idea that for large enough
system sizes the model should exhibit translation invariance. While
translation invariance is broken in the presence of open boundary conditions,
this will still provide a reasonable approximation. In the second phase, each
parameter can be independently varied in a VQE. At this stage, we allow for an
extremely large number of 100,000 iterations to assess the best possible
performance of each ansatz. In practice, we typically observe convergence
after only a few hundred iterations in most cases. As we increment the layer
count from $k$ to $k+1$, the parameters obtained from the prior warm-start
phase for $k$ layers are used again in the first $k$ layers to being the warm
start for $k+1$ layers.

To demonstrate that our VQE ansatz is also suited for current quantum devices,
we take in a second step the parameters obtained from the classical
simulations and perform inference runs on a quantum device. In particular, we
prepare the corresponding states on IBM’s superconducting quantum hardware and
measure the electric field density as well as the particle number. We explore
a wide range of parameter regimes ranging from weak coupling, $m/g\gg 1$, in
which we expect the first-order phase transition to happen up to $m/g\ll 1$,
in which limit we do not expect any transition. In order to obtain the best
possible hardware results, we utilize several error mitigation techniques,
which we outline in the next section.

###  III.2 Error mitigation techniques

In order to mitigate various errors occurring on the quantum devices we use,
we employ a composite mitigation strategy combining several techniques. The
backbone of our composite approach is zero-noise extrapolation (ZNE) [37, 53].
We use digital circuit folding, i.e. inserting pairs of unitaries
$U^{\dagger}U$ in the circuit, which would result in an identity operation on
an ideal quantum computer. For real quantum devices with noise, this
effectively allows for running the same circuit at different (and larger)
noise levels. Subsequently, the results can be extrapolated to zero noise with
an appropriately chosen fitting function [37]. While ZNE allows for addressing
incoherent errors, it cannot mitigate coherent errors.

To counter coherent errors, we employ Pauli twirling to two-qubit CNOT gates
and to measurements, in order to transform their coherent quantum error
channels into incoherent ones [54, 39, 38]. A two-qubit gate can be twirled by
adding random single-qubit Pauli gates before and after it in a way that the
action of the resulting operation is logically unchanged compared to the
original one. Sequential single-qubit gates resulting from twirling
overlapping CNOT gates can be merged to minimize the gate count in the
circuit. Averaging the results over many random instances of twirled circuits
turns coherent errors into incoherent ones. In our case, we use 10 twirls and
average over their result for each noise factor in zero-noise extrapolation.
Measurements can be twirled in a similar manner, except that a random single-
qubit Pauli is inserted only before the measurement. The logical change of
basis is then corrected in post-processing. Different instances of twirled
circuits can be generated by sampling new single-qubit Pauli gates.

Lastly, dynamical decoupling is harnessed to tackle decoherence. During
extended idling periods comparable to the decoherence time, crosstalk between
neighboring qubits can cause decoherence. By leveraging externally
controllable interactions, such as spin flip cycles, during these idle
intervals, we can suppress this form of decoherence. Here, we use the $XX$
sequence, or two X pulses placed uniformly in all idle times [40].

###  III.3 Extrapolation to the continuum limit with Matrix Product States

In order to estimate the system sizes required to obtain a faithful continuum
limit and to examine the behavior of both fermion discretizations towards that
limit, we need to reach system sizes that cannot directly be simulated on a
classical computer. To go beyond the limit of a few tens of qubits, we use
MPS. MPS are an entanglement-based ansatz for quantum states [55, 56, 57, 58]
that can efficiently represent ground states of local gapped Hamiltonians
[59]. For a system of $N$ sites with open boundary conditions, the MPS ansatz
reads

$\displaystyle\ket{\psi}=\sum_{\sigma_{0},..,\sigma_{N-1}}A^{\sigma_{0}}_{0}A^{\sigma_{1}}_{1}\dots
A^{\sigma_{N-1}}_{N-1}\ket{\sigma_{0}}\otimes\ket{\sigma_{1}}...\otimes\ket{\sigma_{N-1}},$
(19)

where $\\{\ket{\sigma_{k}}\\}_{k=1}^{d}$ is a local basis for each site, and
the $A^{\sigma_{k}}_{k}$ are complex $D\times D$ matrices for $0<k<N-1$. The
first (last) tensor $A^{\sigma_{0}}_{0}$ ($A^{\sigma_{N-1}}_{N-1}$) is a
$D$-dimensional row (column) vector. For a fixed combination of physical
indices, the coefficient of the wave function is parameterized as a product of
matrices, hence the name MPS. The size of the matrices $D$ determines the
amount of variational parameters that are present in the ansatz and limits the
maximum entanglement in the state [55, 56, 57, 58].

The MPS approximation for the ground state of a given Hamiltonian can be found
with standard variational algorithms. The tensors are updated iteratively one
by one, while keeping all others fixed. The optimal tensor at each step is
obtained by finding the ground state of an effective Hamiltonian that
describes the interaction of the tensor with its environment. Repeating the
procedure starting from one of the boundaries and sweeping back and forth
until the relative change of the energy is below a certain tolerance, the
resulting MPS represents an approximation for the ground state of the
Hamiltonian.

MPS algorithms routinely deal with several hundreds of spins and allow us to
check the system sizes that are required for obtaining a reliable continuum
extrapolation for both fermion discretizations. In particular, we use systems
with $N\in\\{70,80,90,100\\}$ physical sites and focus on the continuum limit
at a dimensionless lattice volume $N/\sqrt{x}=30$. These system sizes seem to
be within reach on near-term quantum devices.

##  IV Results

Having summarized the methods in Sec. III that we have used to explore the
lattice Schwinger model, we examine the performance of the ansätze for the VQE
using classical simulations. Furthermore, we demonstrate the first-order phase
transition with quantum hardware using the ground states obtained from the
classcial simulation of the VQE. Finally, we perform the extrapolation to the
continuum with data obtained from an MPS calculation and estimate the resource
required for a reliable continuum limit.

###  IV.1 VQE results

In order to explore the capacity of our ansatz of faithfully representing the
ground states, the VQE is simulated classically in the absence of any noise
following the procedure outlined in the previous section. We focus on two
regimes, on the one hand a large value $m_{\text{lat}}/g=10$, for which we
expect the first-order quantum phase transition to occur, and on the other
hand $m_{\text{lat}}/g=0$, which should be considerably below the critical
mass and, thus, no transition occurs as we change $l_{0}$. For the rest of
this section, we consider a fixed lattice volume of $N/\sqrt{x}=30$.

Focusing on the case of case of Wilson fermions first, we generally observe a
similar performance of both the brick and the ladder ansatz. The choice of
SO(4) gates mostly leads to a slightly better fidelity than the $R_{XX+YY}$
gates, independently of the mass. Since the brick ansatz has a smaller CNOT-
depth than the ladder one for the same number of layers, we deem it more
suitable for near-term quantum hardware, and we focus on the brick ansatz with
SO(4) gates for the Wilson case. Using the two product states described in Sec
III.1 as initial states before and after the transition for the Wilson
discretization with $m_{\text{lat}}/g=10$, the VQE was able to achieve
fidelities above 0.99 with a single layer of the brick ansatz with SO(4) gates
for all system sizes we study. Figure 3 shows the results from the VQE for
electric field density, $\langle L_{W}\rangle$, (first row of Fig. 3) and the
particle number, $\langle P_{W}\rangle$, (third row of Fig. 3) as a function
of $l_{0}$ for system sizes $N=3$, $4$, $5$, $6$, which correspond to 6, 8, 10
and 12 qubits for Wilson fermions. In particular, the results for the VQE show
the distinct discontinuity in the electric field density and the particle
number, indicating the remnant of the first-order phase transition. The
location of the transition is shifted to values of $l_{0}$ much larger than
the prediction for the continuum model, $l_{0}=1/2$. As we increase $N$, we
observe that the transition point moves closer towards smaller values of
$l_{0}$. This behavior can be explained with the effects of the finite volume
and lattice spacing on the transition point, which we discuss in detail in
Appendix D. More specifically, since we increase $N$ and keep the volume
constant, we effectively reduce the lattice spacing, hence going closer to the
continuum result. Comparing the results from the VQE to the ones obtained from
exact diagonalization, we observe excellent agreement, reflecting the high
fidelities that were reached. This can be accounted for by the fact that for
$m_{\text{lat}}/g=10$ the true ground states before and after the transition
are close to the product states we use as initial states. As a result, one
layer of the ansatz is sufficient to generate essentially perfect overlap.

Turning to a smaller lattice mass of $m_{\text{lat}}/g=0$, the ground state is
more complicated, and in general our VQE needs more layers to achieve equally
good fidelities as previously. For the most part, two layers of the brick
ansatz with SO(4) gates gave the best performance resulting in fidelities
above 0.99. In a few occasions, the $R_{XX+YY}$ gates showed a slightly better
performance, and we observed a few instances where a single ansatz layer was
already sufficient to reach fidelities above 0.99. Figure 4 summarizes the
results for the electric field density and the particle number for
$m_{\text{lat}}/g=0$. Compared to the larger mass, we now see that both
quantities exhibit a smooth behavior, as one would expect for masses below the
critical one. It can be seen that for a few values of $l_{0}$, the VQE did not
fully converge to the exact diagonalization, especially for larger values of
$N$ close to $l_{0}=1/2$ (see Figs. 3(b), 3(c), 3(d) and 3(j)). In theses
cases the best fidelity obtained is around 0.85. Looking at the electric field
density and the particle number for larger $N$, we observe a rather larger
slope around $l_{0}=1/2$. This indicates that despite choosing a small lattice
mass, we are not too far away from the first-order phase transition, possibly
due to finite-size effects and the additive mass renormalization. Hence, at
$l_{0}=1/2$ there are presumably two states with similar energy making it
harder for the VQE to converge. Nevertheless, our classically simulated VQE
results are in general in good agreement with the data from exact
diagonalization for a wide range of parameters.

Focusing on the case of staggered fermions and using the initial states and
qubit mappings discussed in Sec. III.1, the ansätze with the $R_{XX+YY}$ gates
generally returned higher fidelities compared to the ones with the $SO(4)$
gates. In the latter case, we observed the tendency that the final parameters
obtained by the optimization routine at the end of the VQE were all close to
$0$, resulting in a state formed only from the initial non-parametric part
(which prepares the $\ket{10}^{\otimes N/2}$ state). In contrast, with
$R_{XX+YY}$ gates both of the architectures performed comparably well and were
able to achieve fidelities above $0.99$. Since the brick architecture has
smaller CNOT-depth, we again focus on this one for the staggered
discretization. Considering $m_{\text{lat}}/g=10$ first, the staggered
discretization behaves qualitatively similar to the Wilson fermions, as Fig. 3
reveals. We again observe the characteristic discontinuities in the electric
field density (second row of Fig. 3) and the particle number (fourth row of
Fig. 3), indicating the first-order quantum phase transition. Comparing to the
Wilson fermions, the transition for the staggered discretization occurs at
larger $l_{0}$, despite the fact that for both discretiaztions we match the
qubit number, which effectively means that for the staggered fermions we have
a larger number of physical sites and, since the volume is fixed to
$N/\sqrt{x}=30$, a finer lattice spacing than for the Wilson fermions. This
effect might be caused by a larger additive mass renormalization compared to
the Wilson case. Also for the staggered discretization, we see the same trend
that with increasing system size it is shifted to smaller values of the
$\theta$-angle (for details, see Appendix D). In particular, the data from
exact diagonalization is in excellent agreement with the one from simulated
VQE with one ansatz layer, thus showing that our ansatz circuit together with
the qubit mapping described in Sec. III is appropriate for capturing the
physics of the model in the regime of large masses.

Going to the opposite regime of small lattice mass, $m_{\text{lat}}/g=0$, we
again have to increase the number of layers to two for the staggered case to
have high fidelities with the exact results, due to the more complicated
nature of the ground state in this regime. The results for two layers of the
brick ansatz are shown in Fig. 4. Similar to the Wilson fermions, also for the
staggered discretization, the electric field density (particle number) now
shows a smooth decrease (increase), indicating that the absence of the first-
order phase transition for the chosen lattice mass. Interestingly, for the
staggered discretization we essentially obtain perfect agreement with the data
from exact diagonalization, and there the results from classically simulated
VQE do not show any deviations for the entire parameter range we study.

In summary, our results for the Wilson case suggest that the brick ansatz with
SO(4) gates with an appropriate initial state is generally favorable, and that
the resources required for approximating the ground state to a good level of
accuracy do not show a strong dependence on the system size in the range we
study. For the staggered case, our data imply the brick ansatz with
$R_{XX+YY}$ gates provides good performance using the same initial state for
the entire parameter regime we explore. Moreover, also for the staggered case
there is no strong dependence on the system size in the range we study.

###  IV.2 Inference runs on quantum hardware

The previous section has shown that our VQE ansatz is able to capture the
relevant ground states efficiently with a small number layers. To demonstrate
the approach is feasible on current and near-term quantum hardware, we perform
inference runs on quantum hardware. To this end, we prepare the ansatz circuit
for the parameters obtained at the end of the classical simulation of the VQE
on a quantum device and measure the electric field density as well as the
particle number. We use IBM’s quantum devices ibm_hanoi, ibm_cusco and
ibm_nazca for our inference runs, where we perform $10^{4}$ measurements in
the computational basis. The results for both fermion discretizations after
applying the error mitigation techniques discussed in Sec. III.2 are shown in
Figs. 3 ($m_{\text{lat}}/g=10$) and in 4 ($m_{\text{lat}}/g=0$), alongside
with the data from the simulated ideal VQE and the results from the exact
diagonalization.

In general, the results from the quantum hardware agree well with the noise-
free simulation of the VQE, except at a few points. This discrepancy may be a
result of several factors. We have observed for example that measurements
performed right after calibration of the quantum hardware were closer to the
simulated result, as compared to measurements done when the last calibration
took place several hours before. Furthermore, the case $m_{\text{lat}}/g=0$
has a few points which do not match well with the corresponding noise-free
ones, compared to the larger mass of $m_{\text{lat}}/g=10$. This can be
related to the fact that the circuits for $m_{\text{lat}}/g=0$ were one layer
deeper, which generally leads to larger effects of noise on the results.

Figure 3: Electric field density $\langle L_{W,S}\rangle$ and particle number
$\langle P_{W,S}\rangle$ against $l_{0}$ with data from quantum hardware
(black crosses), as compared to the noiseless expectation values (red pluses)
and exact diagonalization (blue circles) for Wilson and staggered fermions
respectively. The staggered fermions in this case have the same number of
qubits as the Wilson fermions. Hence, while the title of each column specifies
the $N$ for Wilson fermions, for staggered it is taken to be double that
value. The lattice mass for these data is set to $m_{\text{lat}}/g=10$, so
that we are above the second order phase transition of Fig. 1, without having
to account for the mass shift. Therefore, we can observe the first-order phase
transition. Note that the error bars are much smaller than the y-scale and
thus, are not visible. Figure 4: The description for this figure follows Fig.
3, however the lattice mass here is set to $m_{\text{lat}}/g=0$, so that we
are below the second order phase transition and thus observe no first-order
phase transition as expected. Note that the error bars are much smaller than
the y-scale and thus, are not visible. The blue circles represent the exact
diagonalization, red pluses the noiseless simulations and black crosses the
quantum hardware results.

To assess the performance of the error mitigation procedure, we show in Fig. 5
histograms for the absolute error of the unmitigated and mitigated results for
the particle number and the electric field density. The absolute error (shown
on the $x$-axis) corresponds to the absolute value of the difference between
the ones obtained from the simulated VQE, and the (un)mitigated values from
the quantum device. For the electric field density in the Wilson case and the
particle number in the staggered formulation, the histograms in Figs. 5(a) and
5 (d) show a clear improvement. The mitigated values (blue bars) are
concentrated around smaller errors and have a higher probability of coinciding
with the correct result than the unmitigated ones (red translucent bars). For
the electric field density in the staggered case and the particle number in
the Wilson discretization, we observe a different picture. Figures 5(b) and
5(c) show that the probability of measuring the exact result has decreased
compared to the unmitigated results. Nevertheless, the overall distribution
after mitigation shows a smaller width and is concentrated at smaller values
of the absolute error than the original data. This effect can be explained by
the fact that the ZNE is in general not guaranteed to improve the results, and
occasionally it can make the final result slightly worse (examples of ZNE
extrapolations are shown in detail in Appendix B). In general, we observe a
positive effect of the error mitigation, as the expected error is reduced in
all cases.

Figure 5: Histograms for the absolute errors of the mitigated and unmitigated
data for the electric field density for the Wilson (a) and the staggered
discretization (b) as well as the particle number for the Wilson fermions (c)
and the staggered formulation (d). Blue bars represent the absolute error
between the hardware data after ZNE and the simulated VQE, while red bars
represent the absolute error between unmitigated hardware data and the
simulated VQE. Each subplot includes all $N,m_{\text{lat}}/g,l_{0}$.

Having demonstrated, that our ansatz can be reliably executed on current and
near-term quantum hardware for system sizes corresponding 6-12 qubits, we
examine if a near-term device with $\mathcal{O}(100)$ qubits would allow for
reliably studying the continuum limit of the model.

###  IV.3 Continuum extrapolations with MPS

To investigate the system sizes required for a reliable extrapolation to the
continuum limit as well as to compare the performance of both fermion
discretizations towards that limit, we use MPS simulations (for details, c.f.
Appendix C). We use again a fixed lattice volume of $N/\sqrt{x}=30$, where we
restrict $N\in\\{70,80,90,100\\}$. This allows us to study the fixed volume
continuum limit, $ag\to 0$, with resources that should be accessible on
current and near-term quantum hardware. In order to be able to follow a line
of constant renormalized mass, towards the continuum limit, we determine the
mass shift for each value of $N$ and $l_{0}$ following the ideas in Ref. [26].
Here we focus on a small value, $m_{r}/g=0.01$, as for large masses the ground
state of the model is closer to a product state and small masses are in
general more challenging. In particular, we will focus on the electric field
density, as for this observable our continuum extrapolations can be compared
to the prediction from mass perturbation theory [42], which reads

$\frac{\mathcal{F}}{g}=\frac{e^{\gamma}}{\sqrt{\pi}}\Big{(}\frac{m}{g}\Big{)}\sin{\theta}-8.9139\text{
}\frac{e^{2\gamma}}{4\pi}\Big{(}\frac{m}{g}\Big{)}^{2}\sin{(2\theta)},$ (20)

where $\gamma=0.5772156649$ is the Euler-Mascheroni constant.

Figure 6 shows the continuum result for the electric field density obtained
for each discretization in comparison with prediction from mass perturbation
theory. To compare the performance between the staggered and the Wilson
discretization, there are two options. On the one hand, one can choose the
same value of $x=1/(ag)^{2}$ for both, which results in a staggered
formulation that has half of the number of qubits than the Wilson one. On the
other hand, one can choose the same number of qubits for both discretizations,
which leads to double the number of physical sites for the staggered
formulation, which, given the fixed lattice volumes we work with,
simultaneously means a finer lattice spacing than in the Wilson case. For each
case, in the staggered formulation we additionally compare the results
obtained with using the theoretically predicted values of the mass shift in
Eq. (13), and the ones from measuring the mass shift $\text{MS}_{L}$ following
Ref. [26]. Hence, there are in total four combinations of the two above
considerations, which are all shown in Fig. 6.

In general, we observe the most precise results from both Wilson and staggered
fermion when using $\text{MS}_{L}$ and the same $x$ as Wilson. When the
staggered utilized $\text{MS}_{L}$ but retained an equivalent qubit count to
Wilson, the performance dipped, primarily since resources (like bond dimension
$D$) were consistent with the same $x$ scenario, even if the former had an
expanded $N$. Further, the staggered fermions using the $\text{MS}_{t}$ were
not able to match the best performance, since the $\text{MS}_{t}$ does not
take into account the $l_{0}$ dependence to the MS; this effect is exemplified
in the points with higher $l_{0}$ which are further away from perturbation
theory in Fig. 6. Appendix C includes Table 1 which quantifies this deviation
by showing the absolute distance between the points in Fig. 6 and mass
perturbation theory. An $l_{0}$ dependence in the massless Schwinger model
with staggered fermions was also reported in Fig. 7 of [28]. We emphasize that
the $\text{MS}_{t}$ is derived with periodic boundary conditions while we only
consider OBC.

Figure 6: Electric field density $\mathcal{F}/{g}$ against $l_{0}$ for Wilson
fermions and staggered fermions at $m_{\text{r}}/g=0.01$. The following MSE
values are the mean squared error of each of the data as compared to the
continuum mass perturbation theory prediction for the electric field density
$\mathcal{F}/g$ given in Eq. (20). For staggered fermions with $\text{MS}_{t}$
and same $x$ (number of qubits as Wilson), the MSE is $1.041\times 10^{-6}$
($6.16\times 10^{-7}$), and respectively for $\text{MS}_{L}$ $6.143\times
10^{-9}$ ($2.617\times 10^{-8}$). For Wilson fermions the MSE is $7.926\times
10^{-9}$.

Overall, we conclude that, even with these modest resources, the extracted
data closely aligns with the continuum theory, both qualitatively and
quantitatively. However, disregarding the mass shift significantly skews the
extrapolated data in comparison to the theory (for details, see Appendix C).
Hence, if this procedure is to be approached with digital quantum computing,
the mass shift will play an important role in giving reliable continuum limit
results with relatively small system sizes up to $100-200$ qubits.

##  V Conclusion

In this work, we developed VQE ansätze for the lattice Schwinger model
suitable for studying its phase structure on current and near-term quantum
devices in the presence of a topological $\theta$-term, a regime where
conventional Monte Carlo methods suffer from the sign problem. Using two types
of fermion discretizations, Wilson and staggered fermions, we demonstrated
that for both approaches it is possible to find shallow parametric ansatz
circuits that allow for obtaining high fidelities with the exact solution for
a large range of system sizes, masses and $\theta$-values. In particular, we
demonstrated that we can observe the remnants of the expected phase structure
of the first-order phase transition at $\theta=\pi$ occurring in the regime of
large masses. Moreover, preparing the resulting wave function for up to $12$
qubits on IBM’s quantum hardware, we showed that the characteristic
observables revealing the phase structure, the electric field density and the
particle number, can be measured precisely on quantum hardware when using
state-of-the-art error mitigation techniques, notably zero noise
extrapolation, Pauli twirling and dynamical decoupling. In general, our
simulations on classical and quantum hardware suggest that both fermion
discretizations do not show a large difference in performance for the entire
parameter range we study. The ansatz that performed best for the staggered
discretization in our studies is slightly simpler than the one for the Wilson
fermions. Nevertheless, the results from quantum hardware after error
mitigation showed a similar precision for both cases.

Furthermore, we used MPS simulations to estimate the resources for both
fermion discretizations required to obtain the continuum limit. Our results
show that the continuum values for the electric field density and the particle
number can be reliably obtained with $100$ \- $200$ qubits, which is within
reach with current quantum hardware. In addition, our results demonstrate
universality for the two fermion discretizations, as we explicitly observe the
same continuum limit for the electric field density and the particle number
from both approaches. Comparing both discretizations, for the staggered
fermions we observe reduced lattice artifacts towards the continuum limit,
thus allowing us to obtain a continuum limit with the same precision at
slightly larger lattice spacings compared to Wilson fermions. Moreover, for
both cases we observe that taking the additive mass renormalization into
account significantly improves the continuum extrapolation.

Given that for the staggered case our MPS simulations showed better
convergence towards the continuum limit, and our ansatz circuit for the VQE is
slightly simpler, this suggests that the staggered discretization might be
easier to address on current and near-term quantum hardware. Combining our
results with some of the techniques of Ref. [19], it seems possible to reach
larger lattice sizes on the order of 100 qubits, and for the future we plan to
study the continuum extrapolation directly on quantum hardware. Moreover, it
is an interesting question how both fermion discretizations perform in higher
dimensions. While for the Schwinger model we integrated the gauge fields out,
this is no longer possible for two and more spatial dimensions. This might
change the performance of both discretiatzion and the complexity of a VQE
anatz required to reliably capture the ground state.

###### Acknowledgements.

This work is partly funded by the European Union’s Horizon 2020 Research and
Innovation Programme under the Marie Sklodowska-Curie COFUND scheme with grant
agreement no. 101034267. P.N. acknowledges financial support from the Cyprus
Research and Innovation Foundation under projects “Future-proofing Scientific
Applications for the Supercomputers of Tomorrow (FAST)”, contract no.
COMPLEMENTARY/0916/0048, and “Quantum Computing for Lattice Gauge Theories
(QC4LGT)”, contract no. EXCELLENCE/0421/0019. This work is funded by the
European Union’s Horizon Europe Framework Programme (HORIZON) under the ERA
Chair scheme with grant agreement no. 101087126. 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). A. C. is supported in part by the Helmholtz Association —“Innopool
Project Variational Quantum Computer Simulations (VQCS).” The authors also
thank the QC4HEP Working Group for discussions. IBM, the IBM logo, and ibm.com
are trademarks of International Business Machines Corp., registered in many
jurisdictions worldwide. Other product and service names might be trademarks
of IBM or other companies.

##  Appendix A Derivation of the Hamiltonian

In this appendix we derive the Hamiltonian for Wilson fermions, in particular,
we provide the details of the fermion ordering chosen. Subsequently, we show
the Jordan Wingner transformation and the final Hamiltonians in spin
formulation used in the simulations for both Wilson and staggered fermions.

We turn first to Wilson fermions starting from the free Dirac Hamiltonian on
the lattice which is given by

$H=a\sum_{n}\psi_{n}^{\dagger}\gamma^{0}\left(-i\gamma^{1}\partial_{1}+m_{\text{lat}}\right)\psi_{n},$
(21)

with $a$ the lattice spacing, $n$ specifying the lattice site and
$m_{\text{lat}}$ the lattice mass. The derivative $\partial_{1}\psi_{n}$
appearing in the above corresponds to the symmetric finite differences
$(\psi_{n+1}-\psi_{n-1})/2a$. We add the Wilson term [26, 46],

$\displaystyle\begin{split}\Delta
H_{W}&=-r\frac{a^{2}}{2}\sum_{n}\overline{\psi}_{n}\left(\partial_{1}\right)^{2}\psi_{n}\\\
&=-r\frac{a^{2}}{2}\sum_{n}\overline{\psi}_{n}\left(\frac{\psi_{n+1}+\psi_{n-1}-2\psi_{n}}{a^{2}}\right),\end{split}$
(22)

to $H$, where we have used the spatial discrete symmetrized second derivative
for $(\partial_{1})^{2}$ and the Wilson parameter $r$. For the $\gamma$
matrices we make the choice $\gamma^{0}=X$, $\gamma^{1}=iZ$, with $X,Z$, being
the standard Pauli matrices. To ensure invariance under local symmetry
transformations, we have to introduce the link operators $U_{n}$, and its
canonical conjugate the electric field operator $E_{n}=gL_{n}$, with $g$ being
the coupling between fermions and bosons. The resulting gauge invariant
Hamiltonian reads

$\displaystyle\begin{aligned}
H_{W}&=\sum_{n}\biggr{(}-\overline{\psi}_{n}\left(\frac{r+i\gamma^{1}}{2}\right)U_{n}\psi_{n+1}\\\
&+\overline{\psi}_{n}\left(\frac{-r+i\gamma^{1}}{2}\right)U^{\dagger}_{n}\psi_{n-1}\\\
&+(am_{\text{lat}}+r)\overline{\psi}_{n}\psi_{n}+\frac{a}{2}\left(E_{n}+\frac{g\theta}{2\pi}\right)^{2}\biggr{)},\end{aligned}$
(23)

where we explicitly consider the background electric field $g\theta/2\pi$. The
gauge invariant states have to fulfill Gauss’s law, $L_{n}-L_{n-1}=Q_{n}$,
with the charge operator given by $Q_{n}=a\psi_{n}^{\dagger}\psi_{n}-1$. For
open boundary conditions, this expression can be solved iteratively, where we
find $L_{n}=\varepsilon_{0}+\sum_{k=0}^{n}Q_{k}$ with $\varepsilon_{0}$ the
electric field on the left boundary. Similar to the main text, we consider
here $\varepsilon_{0}=0$. Substituting this into the Hamiltonian and applying
the unitary transformation from Ref. [47], we can fully remove the gauge
fields. For numerical calculations, it is convenient to make the fields $\psi$
dimensionless with $\phi_{n,\alpha}=(-1)^{n}\sqrt{a}\psi_{n,\alpha}$, where
the index $\alpha\in\\{1,2\\}$ specifies the spinor component. We also make
the Hamiltonian dimensionless by rescaling it to $\tilde{W}=2H/ag^{2}$.
Contrary to Refs. [46, 26], we additionally apply the transformation
$\phi_{n,\alpha}\to(-1)^{n}(-i)^{\alpha}\phi_{n,\alpha}$ to the fermionic
fields, which results in a completely real Hamiltonian.

In order translate the fermionic fields to qubits, we choose to map the spinor
$\phi_{n,\alpha}$ linearly using the following mapping
$\phi_{n,\alpha}\to\chi_{2n-\lfloor\frac{\alpha}{2}\rfloor+1}$. Now we can
apply a Jordan-Wigner transformation [48],
$\chi_{n}=\prod_{k<n}(iZ_{k})\sigma^{-}_{n}$. The charge operator after this
transformation becomes $Q_{n}=(Z_{2n}+Z_{2n+1})/2$, with $n\in\\{0,N-1\\}$. To
enforce vanishing total charge, in case it is needed, we use a penalty term
$\Lambda$ to the Hamiltonian, $W_{W}=\tilde{W}+\Lambda$, where
$\Lambda=\lambda\left(\sum_{n=0}^{N-1}Q_{n}\right)^{2}$. Provided $\lambda$ is
chosen large enough, we can ensure we remain in the sector with vanishing
total charge. A few algebraic manipulations on nested sums then results in the
final Hamiltonian

$\displaystyle W_{W}$
$\displaystyle=x(r-1)\sum_{n=0}^{N-2}\left(\sigma^{+}_{2n+1}Z_{2n+2}Z_{2n+3}\sigma^{-}_{2n+4}+\text{h.c.}\right)$
$\displaystyle+\frac{x(r+1)}{2}\sum_{n=0}^{N-2}\left(X_{2n+1}X_{2n+2}+Y_{2n+1}Y_{2n+2}\right)$
$\displaystyle+\left(\frac{m_{\text{lat}}}{g}\sqrt{x}+xr\right)\sum_{n=0}^{N-1}\left(X_{2n+1}X_{2n+2}+Y_{2n+1}Y_{2n+2}\right)$
$\displaystyle+\frac{1}{2}\sum_{n=0}^{2N-1}\sum_{k=n+1}^{2N-1}\left(N-\left\lceil\frac{k+1}{2}\right\rceil+\lambda\right)Z_{n}Z_{k}$
(24)
$\displaystyle+l_{0}\sum_{n=0}^{2N-3}\left(N-\left\lceil\frac{n+1}{2}\right\rceil\right)Z_{n}$
$\displaystyle+l_{0}^{2}(N-1)+\frac{1}{4}N(N-1)+\frac{\lambda N}{2},$

where $\lceil\cdot\rceil$ is the ceiling function to the nearest integer
greater than the input.

For the staggered case, we directly apply the Jordan-Wigner transformation to
Eq. (15) and set $\phi_{n}=\prod_{k<n}(iZ_{k})\sigma^{-}_{n}$. Similar to the
Wilson case, we add the penalty term to restrict the Hamiltonian to the
subspace of vanishing total charge. After the Jordan-Wigner transformation,
for staggered fermions, the charge operator becomes $Q_{n}=(Z_{n}+(-1)^{n})/2$
and the final Hamiltonian including the penalty term is given by

$\displaystyle\begin{aligned}
W_{S}&=\frac{x}{2}\sum_{n=0}^{N-2}\left(X_{n}X_{n+1}+Y_{n}Y_{n+1}\right)\\\
&+\frac{1}{2}\sum_{n=0}^{N-2}\sum_{k=n+1}^{N-1}(N-k-1+\lambda)Z_{n}Z_{k}\\\
&+\sum_{n=0}^{N-2}\left(\frac{N}{4}-\frac{1}{2}\left\lceil\frac{n}{2}\right\rceil+l_{0}(N-n-1)\right)Z_{n}\\\
&+\frac{m_{\text{lat}}}{g}\sqrt{x}\sum_{n=0}^{N-1}(-1)^{n}Z_{n}\\\
&+l_{0}^{2}(N-1)+\frac{1}{2}l_{0}N+\frac{1}{8}N^{2}+\frac{\lambda}{4}N.\end{aligned}$
(25)

##  Appendix B Zero noise extrapolations

Here we discuss the technical details of the ZNE that we used for our
inference runs. In our implementation, global folding was applied with noise
factors 1, 3, and 5, with 10,000 shots collected from each circuit. The noise
factors are defined as follows: if the circuit is represented by the unitary
$U$, then at noise factor $f$, where $f$ is an odd number, a total of
$(f-1)/2$ pairs of unitaries $U^{\dagger}U$ is added to the circuit. For
example, noise factor 3, the circuit would be $U(U^{\dagger}U)$ and at noise
factor 5 it would be $U(U^{\dagger}U)(U^{\dagger}U)$.

To get a final zero-noise extrapolated expectation value, we first compute
preliminary zero-noise extrapolated values using three approaches: a linear
fit to noise factors 1, 3, and 5; a linear fit to noise factors 1 and 3; and a
second-order polynomial fit to noise factors 1, 3, and 5. We then calculate
the final zero-noise extrapolated value as a weighted average of the
preliminary values, weighted according to their variance. Their variance
emanates from the fitting procedure performed on points which themselves have
a variance. Each point within a fit, i.e. each point at a given noise factor,
is the average of the 10 twirls used from the Pauli twirling error mitigation
technique described in Sec. III.2. Hence, the corresponding error on each
point within a fit comes from this average. The error of each one of these 10
twirls before twirl averaging is taken to be the square root of the ratio of
the variance of the measurement divided by the number of shots.

In Fig. 7 we show specific parameter examples of the ZNE fits that contribute
to Fig. 5. We present representative cases for both situations, where the ZNE
was able to improve the expectation value of the observable to be measured and
where it was not. Note that Fig. 5 shows that with respect to all experiments
carried out, the ZNE was able to improve the results.

Figure 7: Electric field density $\langle L_{W,S}\rangle$ and particle number
$\langle P_{W,S}\rangle$ against noise factor for zero noise extrapolations
(ZNE) on the data (black crosses) obtained from quantum hardware for different
parameters specified in the titles of each subplot with Wilson and staggered
fermions. The red pluses are noiseless values of the observables and the blue
circles are the exact diagonalization results. For the fits we fit a first-
order polynomial to either all points (blue continuous line) or to points with
noise factor 1, 3 (dashed orange line). We also fit a second order polynomial
to all points (green dotted dashed line). The brown triangle is the weighted
average of the extrapolated points of these 3 fits. The plots on the left
$(a-d)$ are showing examples where the ZNE helped take the expectation values
measured closer to the noiseless result, while on the right $(e-h)$ we show
occasions where the ZNE might not be able to improve the result. However,
overall Fig. 5 indicates that ZNE can significantly improve the results. Where
the error bars are not visible, they are smaller than the markers and y-scale.

##  Appendix C Tensor networks

In order to find an MPS approximation for the ground states of the two
Hamiltonians we consider, we use the two-site variational ground state search
algorithm from the ITensors library [60]. The number of iterations was fixed
to 500 and the tolerance $\eta$ for the relative change in energy between two
sweeps was set to $10^{-9}$, which in turn gives an approximation for the
error on the electric field density given by the value of the electric field
density multiplied by $\sqrt{\eta}$ [61]. Additionally, we restrict the MPS to
the sector of vanishing total charge by applying the symmetry directly to the
ansatz [62], thus no penalty term is used in the spin Hamiltonians.

In order to compare between the Wilson and staggered fermion discretization
approaches, we carry out a continuum limit extrapolation. This implies
following lines of constant physics in the physical phase space, which
requires the calculation of the mass shift for each $N,x,l_{0}$. With our work
confined to a fixed volume of $N/\sqrt{x}=30$, this necessitates the
determination of the mass shift in Eq. (5) for individual values of $N$ and
$l_{0}$ only. In this work, we followd the method based on the electric field
density introduced in Ref. [26]. Table 1 complements Fig. 6 and quantifies the
effect of not accounting for an $l_{0}$ dependence on the mass shift with
staggered fermions. Finally, Fig. 10 is qualitatively showing the effect of
not accounting for the mass shift when using Wilson fermions, which was also
observed in [26]. The bond dimension for all MPS data was set to
$D=20,40,60,80$ and we then extrapolate the last 3 points to infinite $D$ in
the electric field density against $1/D$ plot as shown in Fig. 8. The central
value is taken to be the average of the extrapolated value of the linear fit
and the electric field density value at the highest $D$ calculated, while the
error is approximated to be half the difference between these values. The
smaller $D$ solutions were given as initial ansatz to the higher $D$, while
$D=20$ was initialized with a random MPS.

After calculating the mass shift, we proceed with the calculation of the
electric field density for a fixed renormalized mass of $m_{\text{r}}/g=0.01$
and for different values of $l_{0}$. This procedure enables the comparison of
our continuum extrapolation quantities with the theoretical results from mass
perturbation theory [42]. We use $N$ from 70 to 100 in steps of 10, which in
turn determines the $ag=1/\sqrt{x}$ values through the fixed volume. After
bond dimension extrapolations of electric field density, we extrapolate this
observable to $ag=0$ as shown in Fig. 9. For the Wilson case, we take the
extrapolated value of the linear fit and assign an approximated error to it
which is taken to be the difference between that value and the extrapolated
value of a quadratic fit. For the staggered case, it was observed that second-
order polynomial fits were more suitable and that they converge faster to the
continuum; therefore we keep the extrapolated value of the quadratic fit and
again approximate the error as in the Wilson case. To avoid boundary effects,
the electric field density for Wilson fermions is measured using the middle
link as shown in Eq. (10), and for staggered fermions, the four middle links
are measured to also reduce the staggering effect.

Figure 8: Electric field density $\langle L_{W}\rangle$ against inverse bond
dimension $1/D$ to extrapolate to $D\to\infty$ with $D=40,60,80$. This example
is for Wilson fermions at $N=100$, $l_{0}=0.1$,
$m_{\text{lat}}/g=-0.08236266$. Similar behavior was observed for staggered
fermions as well. Figure 9: Electric field density $\langle L_{W,S}\rangle$
against $ag$ to extrapolate to the continuum limit with fixed volume
$N/\sqrt{x}=30$, fixed $m_{\text{r}}/g=0.01$, using $N=70,80,90,100$. The
continuous orange line is a linear fit and the dashed blue line is a second
order polynomial fit. $(a-c)$ is showing Wilson fermions for
$l_{0}=0.1,0.4,0.6$ respectively and $(d-f)$ the same for staggered fermions
with the same $x$ as Wilson fermions and the electric field density method for
the mass shift [26]. Note that the error bars are much smaller than the
markers and thus, are not visible.

$l_{0}$ | 0.1 | 0.2 | 0.3 | 0.4 | 0.6 | 0.7
---|---|---|---|---|---|---
Wilson | 4.45001467e-05 | 1.00338853e-05 | 5.58001106e-05 | 3.04735487e-05 | 3.64867975e-05 | 2.17138610e-05
| Staggered
---
($\text{MS}_{t}$, same num. qubits as Wilson)}
0.00012783 | 0.00016672 | 0.00022494 | 0.0003207 | 0.00045197 | 0.00046271
| Staggered
---
($\text{MS}_{L}$, same x as Wilson)
6.97181950e-06 | 4.86364083e-06 | 1.14145910e-05 | 4.21571198e-05 | 5.52293937e-05 | 3.33572489e-05
| Staggered
---
($\text{MS}_{t}$, same x as Wilson)
2.91080596e-04 | 2.63562909e-04 | 8.61951831e-05 | 1.03764598e-04 | 6.79235186e-04 | 6.37947680e-04
| Staggered
---
($\text{MS}_{L}$, same num. qubits as Wilson)
2.05654489e-05 | 1.75062880e-05 | 1.54209990e-05 | 1.40152738e-04 | 5.40892433e-05 | 5.13089222e-05

Table 1: This table corresponds to the absolute error with respect to the
results from mass perturbation theory for the data presented in Fig. 6. For
staggered fermions that use the theoretically predicted mass shift,
$\text{MS}_{t}$, the error tends to grow with $l_{0}$ due to the fact that the
$\text{MS}_{t}$ does not account for an $l_{0}$ dependence. Figure 10:
Electric field density $\mathcal{F}/g$ against $l_{0}$ for Wilson fermions
using no mass shift at $m_{\text{r}}/g=m_{\text{lat}}/g=0.01$. The MSE is
calculated as in Fig. 6 and found here to be $7.939\cdot 10^{-4}$.

##  Appendix D Lattice effects on the phase transition

Here we derive an analytical expression for the location of the first-order
phase transition of the lattice Schwinger model focusing on the case of
$m/g>>1$, i.e. the limit where the kinetic term is negligible.

We start by mapping the electric field energy to its lattice version. The
electric field operator can also be made dimensionless using the coupling $g$
between the fermions and gauge fields, i.e. $E\to gL_{n}$. The mapping from a
continuum theory electric field energy to the corresponding lattice version is

$\int_{-\infty}^{\infty}\frac{E^{2}}{2}dx\to
a\sum_{n=0}^{N-2}\frac{g^{2}L_{n}^{2}}{2}.$ (26)

Our derivation will not depend on the specific lattice fermion formulation,
hence the result can be seen as general and applying to both Wilson and
staggered formulations. The key ingredient is the fact at the point of the
first-order phase transition, which we label with $l_{0}^{*}$, the ground
state is two-fold degenerate, as already described in Sec. II.1.

On the left of the phase transition ($l_{0}<l_{0}^{*}$), the ground state has
no charges present and the total electric field is $l_{0}$ on each link.
Hence, the total energy $E_{l_{0}<l_{0}^{*}}$ is just the electric field
energy

$E_{l_{0}<l_{0}^{*}}=a(N-1)g^{2}\frac{l_{0}^{2}}{2},$ (27)

where we have used that on a lattice with open boundaries we have $N-1$ links.
On the right of the phase transition ($l_{0}>l_{0}^{*}$), the ground state has
a negative charge at the left edge of the system and a positive one at the
right edge, connected by an electric flux string lowering total electric field
by one unit. This results in a total electric field of $l_{0}-1$ on each link,
and the total energy consists of the electric field energy plus a mass energy
of $2m_{r}$ from the two charges

$E_{l_{0}>l_{0}^{*}}=a(N-1)g^{2}\frac{(l_{0}-1)^{2}}{2}+2m_{r}.$ (28)

These two energies become equal at $l_{0}=l_{0}^{*}$, and solving for
$l_{0}^{*}$ by equating Eq. (27) with Eq. (28) we find

$l_{0}^{*}=2\frac{m_{\text{r}}}{g}\frac{\sqrt{x}}{N-1}+\frac{1}{2}.$ (29)

In the limit of large $N$, Eq. (29) becomes

$l_{0}^{*}\approx 2\frac{m_{\text{r}}}{g}\frac{\sqrt{x}}{N}+\frac{1}{2},$ (30)

where $\frac{\sqrt{x}}{N}$ is the dimensionless inverse volume. Hence, in the
limit of infinite volume, we obtain the continuum prediction $l_{0}^{*}=1/2$
as expected. If we now take again Eq. (29) and insert for the volume
$N/\sqrt{x}=30$, which we are using in Sec. IV.2 and Sec. IV.3, and express
the renormalized mass in terms of the lattice mass and a mass shift (MS)
according to Eq. (5), we arrive at the result

$l_{0}^{*}=\frac{1}{15}\left(\frac{m_{\text{lat}}}{g}+\text{MS}\right)\frac{1}{1-ag}+\frac{1}{2}.$
(31)

The mass shift decreases with decreasing $ag$ [26, 51]. Further, the factor
$1/(1-ag)$ decreases as $ag\to 0$, which leads to an overall decrease of the
first term and Eq. (31) tends towards the continuum limit. Hence, we see that
the transition point on finite lattices $l_{0}^{*}$ is generally shifted
towars values larger than $1/2$ and gets closer to the continuum theory
prediction for decreasing lattice spacing, as observed in Sec. IV.2 and
supported by Fig. 3.

## References

  * Zhang [2021] H. Zhang, High-throughput design of magnetic materials, Electronic Structure 3, 033001 (2021).
  * Xia  _et al._ [2023] F. Xia, Q. Sun, M. Zhu, Y. Wang, X. Song, Y. Wu, Y. Fang, and W. Li, Formation and magnetic enhancement mechanism of sandwich-structure grain boundary phase in ce magnets, Journal of Magnetism and Magnetic Materials 585, 171112 (2023).
  * Vagov  _et al._ [2023] A. Vagov, T. T. Saraiva, A. A. Shanenko, A. S. Vasenko, J. A. Aguiar, V. S. Stolyarov, and D. Roditchev, Intertype superconductivity in ferromagnetic superconductors, Communications Physics 6, 284 (2023).
  * Court and Cole [2020] C. J. Court and J. M. Cole, Magnetic and superconducting phase diagrams and transition temperatures predicted using text mining and machine learning, npj Computational Materials 6, 18 (2020).
  * Huang and Zhuang [2023] M. Huang and P. Zhuang, Qcd matter and phase transitions under extreme conditions, Symmetry 15, 541 (2023).
  * Callan  _et al._ [1979] C. G. Callan, R. F. Dashen, and D. J. Gross, A theory of hadronic structure, Phys. Rev. D 19, 1826 (1979).
  * Weinberg [1978] S. Weinberg, A new light boson?, Phys. Rev. Lett. 40, 223 (1978).
  * Hang and Wang [2020] X. Hang and J.-P. Wang, Permanent magnet design assisted by antiferromagnet-ferromagnet interface coupling: A monte carlo study, Journal of Magnetism and Magnetic Materials 500, 166360 (2020).
  * Acharyya and Vatansever [2022] M. Acharyya and E. Vatansever, Monte carlo study of the phase diagram of layered xy antiferromagnet, Physica A: Statistical Mechanics and its Applications 605, 128018 (2022).
  * Berthier and Reichman [2023] L. Berthier and D. R. Reichman, Modern computational studies of the glass transition, Nat. Rev. Phys. 5, 102 (2023).
  * mc_[2013] Comparing methods and Monte Carlo algorithms at phase transition regimes: A general overview, in  _Perspectives and Challenges in Statistical Physics and Complex Systems for the Next Decade_, 0 (WORLD SCIENTIFIC, 2013) pp. 3–20.
  * Philipsen [2008] O. Philipsen, Status of Lattice Studies of the QCD Phase Diagram, Progress of Theoretical Physics Supplement 174, 206 (2008).
  * de Forcrand [2009] P. de Forcrand, Simulating QCD at finite density, PoS LAT2009, 010 (2009), arXiv:1005.0539 [hep-lat] .
  * Nagata [2022] K. Nagata, Finite-density lattice QCD and sign problem: Current status and open problems, Prog. Part. Nucl. Phys. 127, 103991 (2022).
  * Funcke  _et al._ [2023] L. Funcke, T. Hartung, K. Jansen, and S. Kühn, Review on Quantum Computing for Lattice Field Theory, PoS LATTICE2022, 228 (2023).
  * Bañuls and Cichy [2020] M. C. Bañuls and K. Cichy, Review on novel methods for lattice gauge theories, Rep. Prog. Phys. 83, 024401 (2020).
  * Bañuls  _et al._ [2020] M. C. Bañuls, R. Blatt, J. Catani, A. Celi, J. I. Cirac, M. Dalmonte, L. Fallani, K. Jansen, M. Lewenstein, S. Montangero, C. A. Muschik, B. Reznik, E. Rico, L. Tagliacozzo, K. V. Acoleyen, F. Verstraete, U.-J. Wiese, M. Wingate, J. Zakrzewski, and P. Zoller, Simulating lattice gauge theories within quantum technologies, The European Physical Journal D 74, 165 (2020).
  * Di Meglio  _et al._ [2023] A. Di Meglio, K. Jansen, I. Tavernelli, C. Alexandrou, S. Arunachalam, C. W. Bauer, K. Borras, S. Carrazza, A. Crippa, V. Croft, _et al._ , Quantum computing for high-energy physics: state of the art and challenges. summary of the qc4hep working group, arXiv:2307.03236 , (2023).
  * Farrell  _et al._ [2023] R. C. Farrell, M. Illa, A. N. Ciavarella, and M. J. Savage, Scalable Circuits for Preparing Ground States on Digital Quantum Computers: The Schwinger Model Vacuum on 100 Qubits, arXiv e-prints , arXiv:2308.04481 (2023), arXiv:2308.04481 [quant-ph] .
  * Martinez  _et al._ [2016] E. A. Martinez, C. A. Muschik, P. Schindler, D. Nigg, A. Erhard, M. Heyl, P. Hauke, M. Dalmonte, T. Monz, P. Zoller, and R. Blatt, Real-time dynamics of lattice gauge theories with a few-qubit quantum computer, Nature 534, 516 (2016).
  * Nguyen  _et al._ [2022] N. H. Nguyen, M. C. Tran, Y. Zhu, A. M. Green, C. H. Alderete, Z. Davoudi, and N. M. Linke, Digital quantum simulation of the schwinger model and symmetry protection with trapped ions, PRX Quantum 3, 020324 (2022).
  * Paulson  _et al._ [2021] D. Paulson, L. Dellantonio, J. F. Haase, A. Celi, A. Kan, A. Jena, C. Kokail, R. van Bijnen, K. Jansen, P. Zoller, and C. A. Muschik, Simulating 2d effects in lattice gauge theories on a quantum computer, PRX Quantum 2, 030334 (2021).
  * Bassman  _et al._ [2021] L. Bassman, M. Urbanek, M. Metcalf, J. Carter, A. F. Kemper, and W. A. de Jong, Simulating quantum materials with digital quantum computers, Quantum Science and Technology 6, 043002 (2021).
  * Thompson and Siopsis [2022] S. Thompson and G. Siopsis, Quantum computation of phase transition in the massive schwinger model, Quantum Science and Technology 7, 035001 (2022).
  * Schwinger [1962] J. Schwinger, Gauge invariance and mass. ii, Phys. Rev. 128, 2425 (1962).
  * Angelides  _et al._ [2023] T. Angelides, L. Funcke, K. Jansen, and S. Kühn, Computing the mass shift of wilson and staggered fermions in the lattice schwinger model with matrix product states, Phys. Rev. D 108, 014516 (2023).
  * Funcke  _et al._ [2020] L. Funcke, K. Jansen, and S. Kühn, Topological vacuum structure of the Schwinger model with matrix product states, Phys. Rev. D 101, 054507 (2020).
  * Azcoiti  _et al._ [2018] V. Azcoiti, E. Follana, E. Royo-Amondarain, G. Di Carlo, and A. Vaquero Avilés-Casco, Massive schwinger model at finite $\theta$, Phys. Rev. D 97, 014507 (2018).
  * Hamer  _et al._ [1982] C. Hamer, J. Kogut, D. Crewther, and M. Mazzolini, The massive Schwinger model on a lattice: Background field, chiral symmetry and the string tension, Nucl. Phys. B 208, 413 (1982).
  * Coleman [1976] S. Coleman, More about the massive schwinger model, Ann. Phys. 101, 239 (1976).
  * Coleman  _et al._ [1975] S. Coleman, R. Jackiw, and L. Susskind, Charge shielding and quark confinement in the massive schwinger model, Ann. Phys. 93, 267 (1975).
  * Honda  _et al._ [2022] M. Honda, E. Itou, Y. Kikuchi, L. Nagano, and T. Okuda, Classically emulated digital quantum simulation for screening and confinement in the schwinger model with a topological term, Phys. Rev. D 105, 014504 (2022).
  * Pederiva  _et al._ [2022] G. Pederiva, A. Bazavov, B. Henke, L. Hostetler, D. Lee, H.-W. Lin, and A. Shindler, Quantum State Preparation for the Schwinger Model, PoS LATTICE2021, 047 (2022).
  * Byrnes  _et al._ [2002] T. M. R. Byrnes, P. Sriganesh, R. J. Bursill, and C. J. Hamer, Density matrix renormalization group approach to the massive schwinger model, Phys. Rev. D 66, 013002 (2002).
  * Peruzzo  _et al._ [2014] 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, Nat. Commun. 5, 4213 (2014).
  * Kitaev [1995] A. Y. Kitaev, Quantum measurements and the Abelian Stabilizer Problem, arXiv e-prints , quant-ph/9511026 (1995), arXiv:quant-ph/9511026 [quant-ph] .
  * Giurgica-Tiron  _et al._ [2020] T. Giurgica-Tiron, Y. Hindy, R. LaRose, A. Mari, and W. J. Zeng, Digital zero noise extrapolation for quantum error mitigation, in  _2020 IEEE International Conference on Quantum Computing and Engineering (QCE)_ (2020) pp. 306–316.
  * van den Berg  _et al._ [2022] E. van den Berg, Z. K. Minev, and K. Temme, Model-free readout-error mitigation for quantum expectation values, Phys. Rev. A 105, 032620 (2022).
  * Cai  _et al._ [2020] Z. Cai, X. Xu, and S. C. Benjamin, Mitigating coherent noise using pauli conjugation, npj Quantum Information 6, 17 (2020).
  * Viola and Lloyd [1998] L. Viola and S. Lloyd, Dynamical suppression of decoherence in two-state quantum systems, Phys. Rev. A 58, 2733 (1998).
  * Buyens  _et al._ [2017a] B. Buyens, S. Montangero, J. Haegeman, F. Verstraete, and K. Van Acoleyen, Finite-representation approximation of lattice gauge theories at the continuum limit with tensor networks, Phys. Rev. D 95, 094509 (2017a).
  * Adam [1997] C. Adam, Massive Schwinger Model within Mass Perturbation Theory, Ann. Phys. 259, 1 (1997).
  * Rothe [2006] H. J. Rothe,  _Lattice gauge theories: an introduction_ (World Scientific Lecture Notes in Physics, 2006).
  * Wilson [1975] K. G. Wilson, Quarks and Strings on a Lattice, in  _13th International School of Subnuclear Physics: New Phenomena in Subnuclear Physics_ (1975).
  * Creutz [1995] M. Creutz, Chiral symmetry on the lattice, Nucl. Phys. B - Proceedings Supplements 42, 56 (1995).
  * Zache  _et al._ [2018] T. V. Zache, F. Hebenstreit, F. Jendrzejewski, M. K. Oberthaler, J. Berges, and P. Hauke, Quantum simulation of lattice gauge theories using Wilson fermions, Quantum Science and Technology 3, 034010 (2018).
  * Hamer  _et al._ [1997] C. J. Hamer, Z. Weihong, and J. Oitmaa, Series expansions for the massive Schwinger model in Hamiltonian lattice theory, Phys. Rev. D 56, 55 (1997).
  * Jordan and Wigner [1928] P. Jordan and E. P. Wigner, About the Pauli exclusion principle, Z. Phys. 47, 631 (1928).
  * Buyens  _et al._ [2017b] B. Buyens, J. Haegeman, F. Hebenstreit, F. Verstraete, and K. Van Acoleyen, Real-time simulation of the schwinger effect with matrix product states, Phys. Rev. D 96, 114501 (2017b).
  * Kogut and Susskind [1975] J. Kogut and L. Susskind, Hamiltonian formulation of wilson’s lattice gauge theories, Phys. Rev. D 11, 395 (1975).
  * Dempsey  _et al._ [2022] R. Dempsey, I. R. Klebanov, S. S. Pufu, and B. Zan, Discrete chiral symmetry and mass shift in the lattice Hamiltonian approach to the Schwinger model, Phys. Rev. Res. 4, 043133 (2022).
  * Liu and Nocedal [1989] D. C. Liu and J. Nocedal, On the limited memory bfgs method for large scale optimization, Mathematical Programming 45, 503 (1989).
  * Majumdar  _et al._ [2023] R. Majumdar, P. Rivero, F. Metz, A. Hasan, and D. S. Wang, Best practices for quantum error mitigation with digital zero-noise extrapolation, arXiv:2307.05203 , (2023).
  * Wallman and Emerson [2016] J. J. Wallman and J. Emerson, Noise tailoring for scalable quantum computation via randomized compiling, Phys. Rev. A 94, 052325 (2016).
  * Orús [2014] R. Orús, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Ann. Phys. 349, 117 (2014).
  * Bridgeman and Chubb [2017] J. C. Bridgeman and C. T. Chubb, Hand-waving and interpretive dance: an introductory course on tensor networks, J. Phys. A: Math. Theor. 50, 223001 (2017).
  * Schollwöck [2011] U. Schollwöck, The density-matrix renormalization group in the age of matrix product states, Ann. Phys. 326, 96 (2011).
  * Bañuls  _et al._ [2020] M. C. Bañuls, R. Blatt, J. Catani, A. Celi, J. I. Cirac, M. Dalmonte, L. Fallani, K. Jansen, M. Lewenstein, S. Montangero, C. A. Muschik, B. Reznik, E. Rico, L. Tagliacozzo, K. Van Acoleyen, F. Verstraete, U.-J. Wiese, M. Wingate, J. Zakrzewski, and P. Zoller, Simulating lattice gauge theories within quantum technologies, The European Physical Journal D 74, 165 (2020).
  * Hastings [2007] M. B. Hastings, An area law for one-dimensional quantum systems, Journal of Statistical Mechanics: Theory and Experiment 2007, P08024 (2007).
  * Fishman  _et al._ [2022] M. Fishman, S. R. White, and E. M. Stoudenmire, The ITensor Software Library for Tensor Network Calculations, SciPost Phys. Codebases , 4 (2022).
  * Haegeman, Jutho [2011] Haegeman, Jutho,  _Variational renormalization group methods for extended quantum systems_ , Ph.D. thesis, Ghent University (2011).
  * Singh  _et al._ [2010] S. Singh, R. N. C. Pfeifer, and G. Vidal, Tensor network decompositions in the presence of a global symmetry, Phys. Rev. A 82, 050301 (2010).