aainstitutetext:  Yukawa Institute for Theoretical Physics, Kyoto University,
Sakyo-ku, Kyoto 606-8502, Japanbbinstitutetext:  Interdisciplinary Theoretical
and Mathematical Sciences Program (iTHEMS), RIKEN,
Wako 351-0198, Japan

# Calculating composite-particle spectra in Hamiltonian formalism and
demonstration in 2-flavor QED${}_{1+1\text{d}}$

Etsuko Itou a,b Akira Matsumoto a and Yuya Tanizaki
itou(at)yukawa.kyoto-u.ac.jp akira.matsumoto(at)yukawa.kyoto-u.ac.jp
yuya.tanizaki(at)yukawa.kyoto-u.ac.jp

###### Abstract

We consider three distinct methods to compute the mass spectrum of gauge
theories in the Hamiltonian formalism: (1) correlation-function scheme, (2)
one-point-function scheme, and (3) dispersion-relation scheme. The first one
corresponds to the conventional Euclidean method in the Monte Carlo
simulations. The second one uses the boundary effect to efficiently compute
the mass spectrum. The third one constructs the excited states and fits their
energy using the dispersion relation with selecting quantum numbers. Each
method has its pros and cons, and we clarify such properties in their
applications to the mass spectrum for the $2$-flavor massive Schwinger model
at $m/g=0.1$ and $\theta=0$ using the density-matrix renormalization group
(DMRG). We note that the multi-flavor Schwinger model at small mass $m$ is a
strongly-coupled field theory even after the bosonizations, and thus it
deserves to perform the first-principle numerical calculations. All these
methods mostly agree and identify the stable particles, pions $\pi_{a}$
($J^{PG}=1^{-+}$), sigma meson $\sigma$ ($J^{PG}=0^{++}$), and eta meson
$\eta$ ($J^{PG}=0^{--}$). In particular, we find that the mass of $\sigma$
meson is lighter than twice the pion mass, and thus $\sigma$ is stable against
the decay process, $\sigma\to\pi\pi$. This is consistent with the analytic
prediction using the WKB approximation, and, remarkably, our numerical results
are so close to the WKB-based formula between the pion and sigma-meson masses,
$M_{\sigma}/M_{\pi}=\sqrt{3}$.

††preprint: YITP-23-98, RIKEN-iTHEMS-Report-23

##  1 Introduction

In recent years, numerical simulations of quantum field theories (QFTs) in the
Hamiltonian formalism have attracted a lot of attention motivated by the rapid
progress of quantum computing technology and also the developments of tensor
network techniques. These methods rely on different disciplines from that of
Monte Carlo simulations for the conventional lattice gauge theories, and thus
they are expected to give complimentary frameworks. One of the remarkable
features is that these methods do not rely on importance sampling, and thus we
may be able to circumvent the issue of sign problems.

With this motivation in mind, we investigate the methods to calculate the mass
spectrum of gauge theories in Hamiltonian formalism. When studying strongly-
coupled QFTs, we often encounter situations where the fundamental degrees of
freedom defining the theory do not appear in the low-energy spectrum. Quantum
chromodynamics (QCD) is a notable, successful example of such phenomena:
Quarks and gluons are confined inside the color-singlet hadrons, and it
explains the physics of strong interaction. The Monte Carlo simulations nicely
predict the hadron spectrum FlavourLatticeAveragingGroupFLAG:2021npn  and also
the physics at finite temperature Borsanyi:2013bia ; HotQCD:2014kol  in the
sign-problem-free regions. Of course, we are currently very far away to
reproduce such tremendous achievements of Monte Carlo simulations, and thus it
is important to develop the counterparts of those calculational techniques in
Hamiltonian formalism.

In this work, we consider three independent methods to compute the mass
spectrum of the $2$-flavor massive Schwinger model using the density-matrix
renormalization group (DMRG):

  * •

correlation-function scheme,

  * •

one-point-function scheme, and

  * •

dispersion-relation scheme.

The first one corresponds to the conventional method using the Euclidean (or
spatial) correlators as in the Monte Carlo simulations. The second one uses
the boundary effect for efficiently computing the mass spectrum, which is
partly motivated by the applications of the Friedel oscillations
PhysRevB.54.13495 ; SHIBATA19971024 . The third one constructs the excited
states and fits their energy using the dispersion relation with selecting
quantum numbers (see, e.g., Refs. Pirvu_2012 ; Haegeman_2012 ; Haegeman_2013
for similar analysis in spin systems). Each method has its pros and cons, and
especially the last one is specific to the Hamiltonian formalism. Our purpose
is to clarify their properties in the concrete studies of the $2$-flavor
massive Schwinger model.

The Schwinger model is a $1+1$d quantum electrodynamics
(QED${}_{1+1\text{d}}$) Schwinger:1962tp , and it is a strongly-coupled theory
like $4$d QCD: The fundamental fermions are confined because of the linear
Coulomb potential, and the low-lying states are composite states like mesons.
Despite its strong-coupling nature, one can calculate many nontrivial aspects
using analytical methods Lowenstein:1971fc ; Casher:1974vf ; Coleman:1975pw ;
Coleman:1976uz ; Manton:1985jm ; Hetrick:1988yg ; Jayewardena:1988td ;
Sachs:1991en ; Adam:1993fc ; Adam:1994by ; Hetrick:1995wq ; Narayanan:2012du ;
Narayanan:2012qf ; Lohmayer:2013eka ; Tanizaki:2016xcu , and this theory has
been used as a benchmark to test new computational methods in previous studies
(see, e.g., Refs. Banuls:2013jaa ; Banuls:2015sta ; Banuls:2016lkq ;
Buyens:2016ecr ; Buyens:2016hhu ; Banuls:2016gid ; Funcke:2019zna ;
Chakraborty:2020uhf ; Honda:2021aum ; Honda:2021ovk ; Honda:2022edn ;
Tomiya:2022chr ; Funcke:2023lli ; Dempsey:2023gib ; Kharzeev:2020kgc ;
deJong:2021wsd ; Nguyen:2021hyk ; Nagano:2023uaq  for numerical studies of
Schwinger model with tensor networks and/or quantum simulations). The mass
spectrum of this model was studied numerically by the Monte Carlo method as
well Fukaya:2003ph , including the region with nonzero $\theta$ angles with
the reweighting technique. As another Monte Carlo based studies, the dual
variable formulations are developed, which successfully eliminates the sign
problem of $1+1$d U(1) gauge theories Gattringer:2015nea ; Gattringer:2015baa
; Gattringer:2018dlw . There is also a numerical approach with the light-cone
quantization in the Hamiltonian formalism Harada:1993va  by using the so-
called Tamm-Dancoff approximation.

The low-energy mass spectrum of the 2-flavor Schwinger model with a theta term
has been studied analytically by using the bosonization technique
Coleman:1976uz ; Hetrick:1995wq . We note, however, that the low-energy
effective theory is strongly coupled if $0<m/g\ll 1$ even after bosonization,
and the details of the prediction rely on some approximations that are not
fully justifiable. It is physically nontrivial if those analytic predictions
are reproduced by the first-principle numerical computations when we go into
the details beyond qualitative aspects.

We performed the DMRG computations using the C++ library of ITensor itensor
to obtain the mass spectrum with the above three methods. We find that all
three methods mostly agree and identify the stable particles, pions $\pi_{a}$
($J^{PG}=1^{-+}$), sigma meson $\sigma$ ($J^{PG}=0^{++}$), and eta meson
$\eta$ ($J^{PG}=0^{--}$), where $J$ denotes the isospin quantum number, $P$
and $G$ is the parity and $G$-parity, respectively. In particular, we observe
that the mass of $\sigma$ meson is lighter than twice the pion mass, and thus
$\sigma$ is stable against the decay process, $\sigma\to\pi\pi$. This implies
that $\sigma$ is a stable particle, not a $\pi\pi$ resonance, and this is a
notable difference compared with $4$d QCD. This is consistent with the
analytic prediction based on the WKB approximation of the Abelian bosonized
description, and our numerical results are so close to the WKB-based formula
between the pion and sigma-meson masses, $M_{\sigma}/M_{\pi}=\sqrt{3}$. Let us
emphasize that this poses an interesting theoretical question on why the
semiclassical prediction works so well even outside its valid regime.

This paper is organized as follows. In Section 2, we review the continuum
2-flavor Schwinger model and the bosonization analysis focusing on the mass
spectrum. In Section 3, we introduce the lattice formulation of the
Hamiltonian and define some observables. In Section 4.1, we briefly explain
our simulation method, the DMRG algorithm. Section 4.2 shows our setup of the
simulation. In Section 5, we present our simulation results for the three
methods. Section 6 is devoted to conclusion and discussion. Appendix A shows
the explicit form of the Hamiltonian and the observables in the spin
representation for DMRG. In Appendix B, we test the validity of the charge
conjugation operator on the lattice in the 1-flavor Schwinger model. In
Appendix C, we discuss the assignment of the flavor index for constructing
MPS. In Appendix D, we investigate how the truncation of the bond dimension
affects the correlation function in the massless 1-flavor Schwinger model.

##  2 Review of the 2-flavor Schwinger model

In this work, we study the 2-flavor Schwinger model, which is
$(1+1)$-dimensional quantum electrodynamics (QED${}_{1+1\text{d}}$) with
$N_{f}=2$ species of Dirac fermion. The Lagrangian density with the Minkowski
metric $\eta_{\mu\nu}=\mathrm{diag}(1,-1)$ is given by

$\mathcal{L}=-\frac{1}{4g^{2}}F_{\mu\nu}F^{\mu\nu}+\frac{\theta}{4\pi}\epsilon_{\mu\nu}F^{\mu\nu}+\sum_{f=1}^{N_{f}}\left[i\bar{\psi}_{f}\gamma^{\mu}\left(\partial_{\mu}+iA_{\mu}\right)\psi_{f}-m\bar{\psi}_{f}\psi_{f}\right],$
(1)

where $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ is the field
strength, $g$ is the gauge coupling, and $\theta$ is the vacuum angle
describing the background electric flux. The index $f$ labels the flavor. We
set the masses of the two fermions equal to $m$.

###  2.1 Global symmetry and composite operators

In the chiral limit ($m=0$), the $2$-flavor Schwinger model has the chiral
symmetry and the $G$-parity symmetry,

$\frac{\mathrm{SU}(2)_{L}\times\mathrm{SU}(2)_{R}}{\mathbb{Z}_{2}}\times(\mathbb{Z}_{2})_{G}\quad(m=0),$
(2)

and the chiral symmetry has an ’t Hooft anomaly. The $G$-parity operation is
the combination of the charge-conjugation with the $\pi$ rotation of the
$\mathrm{SU}(2)_{V}$, which will be discussed later. We note that the
continuous chiral symmetry cannot be spontaneously broken due to the Coleman-
Mermin-Wagner theorem, and the anomaly matching condition is satisfied by the
$\mathrm{SU}(2)$ level-$1$ Wess-Zumino-Witten ($\mathrm{SU}(2)_{1}$ WZW)
conformal field theory. The $\mathrm{SU}(2)_{1}$ WZW model is equivalent to
the self-dual compact boson, and one can explicitly derive it from the
massless $2$-flavor Schwinger model with the Abelian bosonization
Coleman:1976uz .

The massive 2-flavor Schwinger model (1) no longer has the chiral symmetry,
but it maintains the vector-like symmetry,

$\left\\{\begin{array}[]{cc}[\mathrm{SU}(2)_{V}/\mathbb{Z}_{2}]\times(\mathbb{Z}_{2})_{G}&\quad(\theta\in\pi\mathbb{Z}),\\\
\mathrm{SU}(2)_{V}/\mathbb{Z}_{2}&\quad(\text{else}),\end{array}\right.$ (3)

and we call $\mathrm{SU}(2)_{V}/\mathbb{Z}_{2}$ as the isospin symmetry. The
$\mathbb{Z}_{2}$ quotient of $\mathrm{SU}(2)/\mathbb{Z}_{2}$ means that the
local operators always have the integer isospin quantum numbers since the
gauge invariance requires that the local operator must consist of the same
number of $\psi$ and $\bar{\psi}$. We define the isospin operators $J_{a}$ as
conserved charges under this symmetry by

$J_{a}=\frac{1}{2}\int dx\,\bar{\psi}\gamma^{0}\tau_{a}\psi,$ (4)

where $\tau_{a}$ represents Pauli matrices of the isospin space with
$a\in\\{x,y,z\\}$.

When $\theta$ takes some special values, e.g. $\theta=0$, the theory enjoys
the charge conjugation,

$C:A\to-A,\quad\psi\leftrightarrow\mathsf{C}\overline{\psi}^{t},$ (5)

with a suitable element of the Clifford algebra $\mathsf{C}$. We note that
this operation flips the sign of the $\theta$ angle, and thus this symmetry
does not exist for generic values of $\theta$. For general numbers of flavors,
$C$ acts on $\mathrm{SU}(N_{f})/\mathbb{Z}_{N_{f}}$ as an outer automorphism,
i.e. the symmetry group becomes
$[\mathrm{SU}(N_{f})/\mathbb{Z}_{N_{f}}]\rtimes(\mathbb{Z}_{2})_{C}$. When
$N_{f}=2$, however, $\mathrm{SU}(2)$ does not have nontrivial outer
automorphisms, and indeed $C$ just gives the $\pi$ rotation in the isospin
space. Thus, it is convenient to introduce the $G$-parity Gparity ,

$G=Ce^{i\pi J_{y}},$ (6)

so that it commutes with the isospin operation and gives a well-defined
eigenvalue $\pm 1$. Moreover, the $G$ parity acts trivially on the
$\mathrm{SU}(2)_{1}$ WZW theory. Thus, if we find a particle with $G=-1$, we
can immediately tell it remains massive in the chiral limit.

In this paper, we mainly focus on the following composite operators to discuss
the meson spectrum:

$\displaystyle\pi_{a}$ $\displaystyle=-i\bar{\psi}\gamma^{5}\tau_{a}\psi\quad$
$\displaystyle(J^{PG}=1^{-+}),$ (7) $\displaystyle\sigma$
$\displaystyle=\bar{\psi}\psi$ $\displaystyle(J^{PG}=0^{++}),$ (8)
$\displaystyle\eta$ $\displaystyle=-i\bar{\psi}\gamma^{5}\psi$
$\displaystyle(J^{PG}=0^{--}).$ (9)

We call them pion, sigma, and eta operators, respectively, obviously motivated
by the meson spectrum of $4$d QCD. We will often denote $\pi_{3}=\pi$ for
simplicity. Here, we have specified their quantum numbers $J^{PG}$, where $J$
denotes the isospin, and $P$ and $G$ denote the parity and the $G$-parity at
$\theta=0$, respectively. The $(1+1)$d QED is strongly coupled, and it turns
out that the light particles correspond to these operators, and this feature
is reminiscent of $4$d QCD. Here, we would like to note that $\eta$ has
$G=-1$, and thus it remains massive in the chiral limit, which is analogous to
the $U(1)_{A}$ problem Coleman:1975pw ; Coleman:1976uz ; Frohlich:1976mt .
Other mesons, $\pi$ and $\sigma$ have $G=+1$ and actually become massless in
the chiral limit. The massless $\sigma$ particle is an outcome of the absence
of chiral symmetry breaking, unlike the $4$d QCD case.

###  2.2 Phase structure

With $m\not=0$, the system is gapped and has the unique ground state at
generic values of $\theta$. In $(1+1)$d, there is no stable topologically-
ordered state, and the unique gapped ground states are then classified as the
symmetry-protected topological (SPT) states PhysRevB.83.035107 ;
Kapustin:2014gma ; Kapustin:2014tfa . This perspective provides a very
powerful tool to understand the phase structure of the $2$-flavor Schwinger
model.

Let us recall that the massive Schwinger model always has the isospin
symmetry, $\mathrm{SU}(2)_{V}/\mathbb{Z}_{2}$. We can then calculate the
partition function under the presence of the background gauge fields for the
isospin symmetry. Compared with the $\mathrm{SU}(2)$ gauge field, the
$\mathrm{SU}(2)/\mathbb{Z}_{2}$ gauge field has milder cocycle conditions,
which is controlled by the $\mathbb{Z}_{2}$ $2$-form gauge field $w_{2}$ in
addition to the familiar $1$-form gauge field. As a result, at generic values
of $\theta$, the partition function with the background gauge field is
described by the low-energy effective topological action,

$\mathcal{Z}_{\theta}\simeq\exp(i\pi k\int w_{2}),$ (10)

with some $k\sim k+2$. We note that $k$ is a discrete label that distinguishes
the SPT states protected by $\mathrm{SU}(2)_{V}/\mathbb{Z}_{2}$, and thus it
cannot be changed under the continuous change of coupling constants unless
quantum phase transitions happen.

We can prove that two partition functions at $\theta$ and $\theta+2\pi$ are
different in the presence of the background gauge fields. The anomalous
relation can be summarized as Misumi:2019dwq

$\mathcal{Z}_{\theta+2\pi}=\exp\left(i\pi\int
w_{2}\right)\mathcal{Z}_{\theta}.$ (11)

The label $k$ is changed as $k\mapsto k+1$ as we change the $\theta$ angle
from $\theta$ to $\theta+2\pi$, and there must be a phase transitions
separating the $k=0,1$ ground states. It is somewhat customary to assign the
$k=0$ state for $-\pi<\theta<\pi$ and the $k=1$ state for $\pi<\theta<3\pi$,
while, precisely, this assignment depends on the UV-regularization scheme. We
note that the whole story here is quite parallel to that of the anti-
ferromagnetic Heisenberg chain or the $(1+1)$d $\mathbb{C}P^{1}$ sigma model
Haldane:1983ru ; Affleck:1986pq ; Haldane:1988zz ; Affleck:1987vf ;
Komargodski:2017dmc ; Komargodski:2017smk ; Lajko:2017wif ; Tanizaki:2018xto
(except for the properties at $\theta=\pi$ Coleman:1976uz ; Dempsey:2023gib ).

The distinction between the states at $\theta=0$ and $\theta=2\pi$ becomes
more vivid when we take the open boundary condition. Turning on non-zero
$\theta$ corresponds to introducing a background electric field with a
constant magnitude $\theta/2\pi$. When we increase $\theta$ beyond $\pi$, the
background field becomes larger than $1/2$. Then the Dirac fermions with
charges $\pm 1$ are excited at the boundaries to cancel the background field
as much as possible. As a consequence, these boundary states have isospin
$1/2$, which is the projective representation of
$\mathrm{SU}(2)_{V}/\mathbb{Z}_{2}$. This is nothing but the signature of the
nontrivial SPT state. If the system size is large enough, the interaction
between the boundary states is exponentially suppressed, so the two
independent degrees of freedom with isospin $1/2$ yield $2\times 2$ degeneracy
of the ground state. We will use this boundary excitation as a source of iso-
triplet particles to determine the mass of the pion from the one-point
function in Section 5.2. Note that if we increase $\theta$ further beyond
$3\pi$, the boundary excitations become bound states of Dirac fermions with
isospin 1, which can be completely screened by the gauge-invariant particles
inside the bulk, and the ground state would be unique again with the open
boundary condition.

###  2.3 Mass spectrum

In this subsection, we are going to give a relatively detailed review on
analytic predictions about the mass spectrum. There is already a huge effort
for the analytic studies of the multi-flavor Schwinger model, and thus the
reader may wonder how one could obtain something new with numerical studies.
We would like to clarify what kinds of approximations are used in the previous
studies and give the justification of this work on physical aspects.

There are two exactly-solvable limits of the multi-flavor Schwinger model:

  * •

In the heavy fermion limit $m\to\infty$, the model becomes the pure $U(1)$
gauge theory.

  * •

In the chiral limit $m=0$, the model becomes the $SU(N_{f})_{1}$ WZW conformal
field theory and one massive free boson, and they are completely decoupled.

It is then natural to consider perturbations from these limits in order to
investigate the cases of general mass $m$. When $m\gg g$, one can perform
systematic perturbations to study the spectrum. In the opposite case $0<m\ll
g$, however, the systematic perturbation works only for the $1$-flavor case,
and further approximations are necessary for $N_{f}\geq 2$.

By applying the Abelian bosonization for $N_{f}=2$, the fermions are mapped to
the $2\pi$-periodic scalar fields $\phi_{1},\phi_{2}$. The Lagrangian (1) is
then completely equivalent to

$\displaystyle\mathcal{L}$
$\displaystyle=\frac{1}{2g^{2}}F_{01}^{2}+\frac{1}{2\pi}(\phi_{1}+\phi_{2}+\theta)F_{01}$
$\displaystyle+\frac{1}{8\pi}\left((\partial\phi_{1})^{2}+(\partial\phi_{2})^{2}\right)+Cm\rho
N_{\rho}[\cos(\phi_{1})+\cos(\phi_{2})],$ (12)

where $C=e^{\gamma}/(2\pi)$ is a numerical constant, and $N_{\rho}$ denotes
the normal ordering for the contraction with a free field propagator of mass
$\rho$ Coleman:1974bu .111All the UV divergences from loop diagrams are
removed by this prescription, and the theory is independent of the choice of
$\rho$. For the free theory with mass $\rho$, $N_{\rho}[\bullet]$ becomes the
ordinary normal ordering. We can rigorously integrate out the gauge fields as
the Lagrangian is quadratic in terms of $F_{01}$. Changing the basis of bosons
as $\phi_{1,2}=\sqrt{2\pi}\eta-\frac{\theta}{2}\pm\varphi$, the effective
Lagrangian becomes

$\displaystyle\mathcal{L}_{\mathrm{eff}}[\eta,\varphi]=\frac{1}{2}\left[(\partial\eta)^{2}-\mu^{2}\eta^{2}\right]+\frac{1}{4\pi}(\partial\varphi)^{2}+2Cm\rho
N_{\rho}\left[\cos\left(\sqrt{2\pi}\eta-\frac{\theta}{2}\right)\cos(\varphi)\right],$
(13)

where $\mu^{2}=2g^{2}/\pi$. We have the $\mathbb{Z}_{2}$ symmetry acting only
on $\eta$ when $\theta=0$, and this is the $G$-parity. When $m=0$, the massive
$\eta$ and the massless $\varphi$ decouple, as advocated above.

Let us now turn on the small mass, $0<m\ll g$. The $\eta$ particle has the
mass $\mu+O(m)$, but it is hard to compute the $O(m)$ correction due to the
potential infrared divergence in the loop diagrams with the $\varphi$ fields
Coleman:1976uz . Instead, we integrate out $\eta$ at the tree level to discuss
the physics of $\pi$ and $\sigma$ mesons, which gives $\left\langle
N_{\rho}\left[\cos\left(\sqrt{2\pi}\eta-\frac{\theta}{2}\right)\right]\right\rangle=\sqrt{\frac{\mu}{\rho}}\cos\frac{\theta}{2}$
for the free massive $\eta$ with mass $\mu$. The effective theory for
$\varphi$ becomes the sine-Gordon model,

$\mathcal{L}_{\mathrm{SG}}[\varphi]=\frac{1}{4\pi}(\partial\varphi)^{2}+2Cm\cos\frac{\theta}{2}(\mu\rho)^{1/2}N_{\rho}[\cos\varphi].$
(14)

The isospin $\mathrm{SU}(2)_{V}/\mathbb{Z}_{2}$ symmetry is not manifest at
all in the Lagrangian, but it secretly exists quantum mechanically. In
particular, the $z$-component isospin current is given by
$J_{z}^{\mu}=\frac{1}{2\pi}\varepsilon^{\mu\nu}\partial_{\nu}\varphi$, and its
charge $J_{z}=\int dx\frac{1}{2\pi}\partial_{x}\varphi$ counts the winding
number of $\varphi\sim\varphi+2\pi$.

We would like to emphasize that the effective theory (14) is strongly coupled,
and we cannot solve it in the ordinary perturbation for small but nonzero $m$.
What is actually done in the previous literature is the optimized
perturbation; we optimize the renormalization scale $\rho$ so that the
coefficient of the $\cos(\varphi)$ potential becomes $O(\rho^{2})$, and we get

$\rho_{\mathrm{optimzied}}\sim\left|m\sqrt{\mu}\cos(\theta/2)\right|^{2/3}.$
(15)

This is identified as the mass gap caused by the mass perturbation, and this
formula gives the $\theta$-dependence of the lightest meson mass, i.e.
$M_{\pi}$.

The spectrum of the sine-Gordon model was studied by using WKB approximation
Dashen:1975hd . Introducing an extra parameter controlling the kinetic term as
$\frac{1}{4\pi\beta^{2}}(\partial\varphi)^{2}$, the quantum scaling dimension
of $\cos\varphi$ becomes $\Delta=\beta^{2}/2$, so the semiclassical
approximation is valid when $\beta^{2}\to 0$. The model has the soliton and
antisoliton, and let us denote their mass as $M_{\mathrm{SG}}$. Then, Dashen
et al. Dashen:1975hd  predicted the masses of soliton-antisoliton bound states
as

$M_{\mathrm{SG}}^{(n)}=2M_{\mathrm{SG}}\sin\left(\frac{\pi}{2}\frac{n}{4/\beta^{2}-1}\right),$
(16)

with $n=1,2,\cdots<(4/\beta^{2}-1)$. Even though it is subtle if the WKB works
at the self-dual point $\beta^{2}=1$, Coleman got an intriguing observation
using this semiclassical formula Coleman:1976uz .

The nontrivial check for its validity is the recovery of the
$\mathrm{SU}(2)_{V}/\mathbb{Z}_{2}$ symmetry at $\beta^{2}=1$. Substituting
$n=1$ with $\beta^{2}=1$ into (16), we get

$M_{\mathrm{SG}}^{(1)}=2M_{\mathrm{SG}}\sin\frac{\pi}{6}=M_{\mathrm{SG}}.$
(17)

This shows that the lightest soliton-antisoliton bound state has the same mass
as the soliton or antisoliton itself. The soliton and antisoliton have
$J_{z}=\pm 1$, and the soliton-antisoliton bound state has $J_{z}=0$. The
$G$-parity does not act on $\varphi$, so all these states have $G=+1$. Thus,
these three states form the isospin triplet $J^{PG}=1^{-+}$ of mass
$M_{\mathrm{SG}}$, which is identified as the pion in the Schwinger model, and
then $M_{\mathrm{SG}}=M_{\pi}\sim(m\sqrt{\mu}\cos(\theta/2))^{2/3}$.

The mass of the second soliton-antisoliton bound state is given by

$M_{\mathrm{SG}}^{(2)}=2M_{\mathrm{SG}}\sin\frac{\pi}{3}=\sqrt{3}M_{\mathrm{SG}}.$
(18)

This state has $J_{z}=0$, and there is no other state with the same mass. We
then identify it as the $\sigma$ meson in the Schwinger model with
$J^{PG}=0^{++}$. Thus, the semiclassical method predicts that the masses of
pion and sigma meson satisfy

$\displaystyle M_{\sigma}=\sqrt{3}M_{\pi}.$ (19)

Importantly, $M_{\sigma}<2M_{\pi}$. Unlike the $4$d QCD, $\sigma$ is a stable
particle, not a resonance, because the decay $\sigma\to\pi\pi$ is
energetically prohibited.

As we discussed above, the low-energy mass spectra can be estimated by
bosonization. However, it relies on the optimized perturbation and also on the
semiclassical method, and these analyses are not necessarily fully justified.
It is still difficult to compute the exact $m$-dependence or to find other
states with higher energies than $\mu$. Thus, it is worth studying the mass
spectrum of the 2-flavor model by first-principle numerical methods.

##  3 Lattice formulation of the 2-flavor Schwinger model

In this section, we explain the Hamiltonian formalism of the 2-flavor
Schwinger model and its lattice regularization as a generalization to
$N_{f}=2$ from previous research Chakraborty:2020uhf ; Honda:2021aum ;
Honda:2021ovk ; Honda:2022edn . We also define various local and global
observables used in the analysis.

###  3.1 Hamiltonian

First, we introduce the continuum Hamiltonian of the $N_{f}$-flavor Schwinger
model. By introducing a conjugate momentum
$\Pi=\frac{1}{g^{2}}\partial_{0}A^{1}+\frac{\theta}{2\pi}$, the Hamiltonian is
given by

$H=\int
dx\,\left\\{\frac{g^{2}}{2}\left(\Pi-\frac{\theta}{2\pi}\right)^{2}+\sum_{f=1}^{N_{f}}\left[-i\bar{\psi}_{f}\gamma^{1}\left(\partial_{1}+iA_{1}\right)\psi_{f}+m\bar{\psi}_{f}\psi_{f}\right]\right\\}.$
(20)

In Hamiltonian formalism, the physical Hilbert space is constrained by the
Gauss law condition

$\partial_{1}\Pi+\sum_{f=1}^{N_{f}}\psi_{f}^{\dagger}\psi_{f}=0.$ (21)

The electric field corresponds to $E:=\dot{A}_{1}=g^{2}(\Pi-\theta/2\pi)$.
Thus, the theta angle $\theta$ plays the role of the background electric
field. In the periodic boundary condition, the Hamiltonians at $\theta$ and
$\theta+2\pi$ are unitary equivalent, $H_{\theta+2\pi}=U^{\dagger}H_{\theta}U$
with $U=\exp(-i\int A_{1}dx)$, which realizes the $2\pi$ periodicity of
$\theta$.

Next, we consider the lattice regularization of the Hamiltonian. Here we
employ the staggered fermion to define fermions on the lattice Kogut:1974ag ;
Susskind:1976jm . The staggered fermions $\chi_{f,n}$ with the lattice spacing
$a$ represents the discretization of the two-component Dirac fermions
$\psi_{f}(x)$ with the lattice spacing $2a$. The single-component fermions
$\chi_{f,n}$ at the site $n=0,1,\cdots,N-1$ correspond to each component222
The labels $u$ and $d$ of $\psi_{f}$ denote the upper and lower spinor
component respectively. They are nothing to do with up and down quark in QCD
here. of $\psi_{f}(x)$ depending on $n$ as

$\psi_{f}(x)=\begin{pmatrix}\psi_{u,f}(x)\\\
\psi_{d,f}(x)\end{pmatrix}\leftrightarrow\frac{1}{\sqrt{2a}}\begin{pmatrix}\chi_{f,2[n/2]}\\\
\chi_{f,2[n/2]+1}\end{pmatrix}.$ (22)

The number of staggered fermions for each flavor is equal to $N$, thus at each
site, there are two staggered fermions. In this work, we set $N$ to be an even
number.

The gauge field is encoded to U(1) variables $U_{n}\sim\exp(-iaA^{1}(x))$,
defined on the link between the $n$-th and $(n+1)$-th sites, and the conjugate
momentum is replaced by $L_{n}\sim-\Pi(x)$, defined on the $n$-th site. The
canonical commutation relations are given by

$\\{\chi_{f,n}^{\dagger},\,\chi_{f^{\prime},m}\\}=\delta_{ff^{\prime}}\delta_{nm},$
(23)
$\\{\chi_{f,n},\,\chi_{f^{\prime},m}\\}=\\{\chi_{f,n}^{\dagger},\,\chi_{f^{\prime},m}^{\dagger}\\}=0,$
(24) $[U_{n},\,L_{m}]=\delta_{nm}U_{n}.$ (25)

Note that the roles of the staggered fermion operators depend on the site $n$:

$\chi_{f,n}^{\dagger}:\begin{cases}\textrm{creation op. of
particle}&n:\textrm{even}\\\ \textrm{annihilation op. of anti-
particle}&n:\textrm{odd}\end{cases},$ (26)
$\chi_{f,n}:\begin{cases}\textrm{annihilation op. of
particle}&n:\textrm{even}\\\ \textrm{creation op. of anti-
particle}&n:\textrm{odd}\end{cases}.$ (27)

Thus, the operator $\chi_{f,n}^{\dagger}\chi_{f,n}$ counts the number of
particles on the even sites, whereas $\chi_{f,n}\chi_{f,n}^{\dagger}$ counts
the number of anti-particles on the odd sites. Considering that the particle
has an electric charge of $+1$ and the anti-particle has $-1$, the charge
density operator at the site $n$ is given by

$\rho_{f,n}=\chi_{f,n}^{\dagger}\chi_{f,n}+\frac{(-1)^{n}-1}{2}=\begin{cases}\chi_{f,n}^{\dagger}\chi_{f,n}&n:\textrm{even}\\\
-\chi_{f,n}\chi_{f,n}^{\dagger}&n:\textrm{odd}\end{cases}.$ (28)

In this work, we choose the open boundary condition in order to eliminate the
bosonic degrees of freedom having an infinite dimensional Hilbert space. The
Gauss law (21) is also discretize as

$L_{n}-L_{n-1}=\sum_{f=1}^{N_{f}}\rho_{f,n},$ (29)

where the left-hand side corresponds to the divergence of the electric field
and the right-hand side is the charge density (28).

We set the explicit form of the (1+1)d gamma matrices $\gamma^{0}=\sigma^{3}$,
$\gamma^{1}=i\sigma^{2}$ and $\gamma^{5}=\gamma^{0}\gamma^{1}=\sigma^{1}$.
Using the operators $\chi_{f,n}$, $U_{n}$, and $L_{n}$ introduced above, the
lattice Hamiltonian is given by Funcke:2023lli ; Dempsey:2023gib

$\displaystyle H$
$\displaystyle=J\sum_{n=0}^{N-2}\left(L_{n}+\frac{\theta}{2\pi}\right)^{2}$
$\displaystyle+\sum_{f=1}^{N_{f}}\left[-iw\sum_{n=0}^{N-2}\left(\chi_{f,n}^{\dagger}U_{n}\chi_{f,n+1}-\chi_{f,n+1}^{\dagger}U_{n}^{\dagger}\chi_{f,n}\right)+m_{\mathrm{lat}}\sum_{n=0}^{N-1}(-1)^{n}\chi_{f,n}^{\dagger}\chi_{f,n}\right],$
(30)

where $J=g^{2}a/2$ and $w=1/2a$. Here we replace the mass $m$ of the continuum
theory by

$m_{\mathrm{lat}}:=m-\frac{N_{f}g^{2}a}{8}$ (31)

in the lattice Hamiltonian, following the recent proposal Dempsey:2022nys  for
eliminating $O(a)$ correction. In the continuum theory, the chiral limit $m=0$
has the continuous chiral symmetry, and it contains
$[\mathrm{SU}(N_{f})_{V}/\mathbb{Z}_{N_{f}}]\times(\mathbb{Z}_{N_{f}})_{L}$ as
a subgroup. With the above replacement, the lattice theory at $m=0$ maintains
the discrete chiral symmetry $\mathbb{Z}_{2}\subset(\mathbb{Z}_{N_{f}})_{L}$
for even $N_{f}$, and this is the point protected by the remnant of the chiral
symmetry.

By adding up the lattice Gauss law equation (29) from the boundary $n=0$ to
the site $n$, we find that $L_{n}$ can be replaced by

$\displaystyle L_{n}$
$\displaystyle=L_{-1}+\sum_{f=1}^{N_{f}}\sum_{k=0}^{n}\rho_{f,k}$
$\displaystyle=\sum_{f=1}^{N_{f}}\sum_{k=0}^{n}\chi_{f,k}^{\dagger}\chi_{f,k}+\frac{N_{f}}{2}\left(\frac{(-1)^{n}-1}{2}-n\right),$
(32)

where we set $L_{-1}=0$. Furthermore, we can set $U_{n}=1$ since the degrees
of freedom of $U_{n}$ can be absorbed by the U(1) phase of $\chi_{n}$. Then
the lattice Hamiltonian is written only by the fermions as
$H=H_{J}+H_{w}+H_{m}$, where the gauge part $H_{J}$ is given by

$H_{J}=J\sum_{n=0}^{N-2}\left[\sum_{f=1}^{N_{f}}\sum_{k=0}^{n}\chi_{f,k}^{\dagger}\chi_{f,k}+\frac{N_{f}}{2}\left(\frac{(-1)^{n}-1}{2}-n\right)+\frac{\theta}{2\pi}\right]^{2},$
(33)

and the kinetic term $H_{w}$ and the mass term $H_{m}$ of the fermions are

$H_{w}=-iw\sum_{f=1}^{N_{f}}\sum_{n=0}^{N-2}\left(\chi_{f,n}^{\dagger}\chi_{f,n+1}-\chi_{f,n+1}^{\dagger}\chi_{f,n}\right),$
(34)
$H_{m}=m_{\mathrm{lat}}\sum_{f=1}^{N_{f}}\sum_{n=0}^{N-1}(-1)^{n}\chi_{f,n}^{\dagger}\chi_{f,n}.$
(35)

###  3.2 Map to the spin system

Now, we map the Hamiltonian written by the staggered fermions to the spin
Hamiltonian. Such a spin Hamiltonian formalism is useful to apply tensor
network methods and quantum computations.

The $N_{f}\times N$ degrees of freedom of the staggered fermion $\chi_{f,n}$
can be described by the same number of spin-1/2 degrees of freedom. The
Hilbert space of such a spin system is given by

$\mathcal{H}=\bigotimes_{f=1}^{N_{f}}\bigotimes_{n=0}^{N-1}\mathcal{H}_{f,n},$
(36)

where $\mathcal{H}_{f,n}$ is the local Hilbert space of the single spin-1/2
state. A general state $\ket{\Psi}$ in this Hilbert space can be described by
a superposition of all possible spin configurations $\bm{s}$,

$\ket{\Psi}=\sum_{\bm{s}}\Psi(\bm{s})\ket{\bm{s}},$ (37)
$\ket{\bm{s}}\in\left\\{\left.\bigotimes_{f=1}^{N_{f}}\bigotimes_{n=0}^{N-1}\ket{s_{f,n}}_{f,n}\right|\ket{s_{f,n}}_{f,n}=\ket{\uparrow},\ket{\downarrow}\right\\}.$
(38)

The spin up $\ket{\uparrow}$ and down $\ket{\downarrow}$ state are the
eigenstates of the Pauli matrix $\sigma^{z}$ with the eigenvalues $+1$ and
$-1$, respectively.

The map to the spin system can be achieved by the so-called Jordan-Wigner
transformation. The fermion operators $\chi_{f,n}$ for the two flavors $f=1,2$
are represented by spin operators as follows:

$\displaystyle\chi_{1,n}$
$\displaystyle=\sigma_{1,n}^{-}\prod_{j=0}^{n-1}(-\sigma_{2,j}^{z}\sigma_{1,j}^{z}),$
$\displaystyle\chi_{1,n}^{\dagger}$
$\displaystyle=\sigma_{1,n}^{+}\prod_{j=0}^{n-1}(-\sigma_{2,j}^{z}\sigma_{1,j}^{z}),$
(39) $\displaystyle\chi_{2,n}$
$\displaystyle=\sigma_{2,n}^{-}(-i\sigma_{1,n}^{z})\prod_{j=0}^{n-1}(-\sigma_{2,j}^{z}\sigma_{1,j}^{z}),$
$\displaystyle\chi_{2,n}^{\dagger}$
$\displaystyle=\sigma_{2,n}^{+}(i\sigma_{1,n}^{z})\prod_{j=0}^{n-1}(-\sigma_{2,j}^{z}\sigma_{1,j}^{z}),$
(40)

where we define

$\sigma_{f,n}^{\pm}=\frac{1}{2}(\sigma_{f,n}^{x}\pm i\sigma_{f,n}^{y}).$ (41)

The Pauli matrices $\sigma_{f,n}^{a}$ ($a=x,y,z$) act on the spin
$\ket{s_{f,n}}_{f,n}$ at the site $n$ of the flavor $f$. They do not commute
only if they are on the same site of the same flavor, so that

$\left[\sigma_{f,n}^{a},\,\sigma_{f^{\prime},n^{\prime}}^{b}\right]=2i\delta_{ff^{\prime}}\delta_{nn^{\prime}}\epsilon^{abc}\sigma_{f,n}^{c}.$
(42)

We can check that the canonical anti-commutation relations (23) and (24) are
satisfied thanks to the properties of the Pauli matrices.

Note that this is not a unique way of translation to the spin system which
realizes the anti-commutation relations. Different transformations give
different representations of the original Hamiltonian. We choose this
transformation since various local operators can be constructed by only a few
numbers of the Pauli matrices. The spin representation of the Hamiltonian and
the observables defined above are summarized in Appendix A.

###  3.3 Local observables

Let us consider the meson operators (7), (9), and (8) on the lattice. Based on
the continuum descriptions, it is natural to define the lattice version of
these operators as follows:

$\pi(n):=PS_{1,n}-PS_{2,n},$ (43) $\eta(n):=PS_{1,n}+PS_{2,n},$ (44)
$\sigma(n):=S_{1,n}+S_{2,n}.$ (45)

Here $S_{f,n}$ and $PS_{f,n}$ are the scalar and pseudo-scalar operators for
the flavor $f=1,2$ on the lattice, respectively. In order to obtain their
explicit form, we rewrite the scalar condensate $(\bar{\psi}\psi)_{f}$ by the
staggered fermion, so that

$\displaystyle(\bar{\psi}\psi)_{f}$
$\displaystyle=(\psi_{u}^{\dagger}\psi_{u}-\psi_{d}^{\dagger}\psi_{d})_{f},$
$\displaystyle=\frac{1}{2a}(-1)^{n}(\chi_{f,n}^{\dagger}\chi_{f,n}-\chi_{f,n+1}^{\dagger}\chi_{f,n+1}).$
(46)

Similarly, the pseudo-scalar condensate $(\bar{\psi}\gamma^{5}\psi)_{f}$ is
given by

$\displaystyle(\bar{\psi}\gamma^{5}\psi)_{f}$
$\displaystyle=(\psi_{u}^{\dagger}\psi_{d}-\psi_{d}^{\dagger}\psi_{u})_{f},$
$\displaystyle=\frac{1}{2a}(-1)^{n}(\chi_{f,n}^{\dagger}\chi_{f,n+1}-\chi_{f,n+1}^{\dagger}\chi_{f,n}).$
(47)

These operators have a site-by-site fluctuation due to the staggered fermion.
Here we define the lattice scalar condensate operator by the two-site average
of (46), namely

$\displaystyle S_{f}(n)$
$\displaystyle:=\frac{1}{2}\left[(\bar{\psi}\psi)_{f,n-1}+(\bar{\psi}\psi)_{f,n}\right],$
$\displaystyle=\frac{1}{4a}(-1)^{n}(-\chi_{f,n-1}^{\dagger}\chi_{f,n-1}+2\chi_{f,n}^{\dagger}\chi_{f,n}-\chi_{f,n+1}^{\dagger}\chi_{f,n+1}),$
(48)

for $n=1,2,\cdots,N-2$. The lattice pseudo-scalar condensate operator is also
defined by the two-site average of (47) with a factor $-i$,

$\displaystyle PS_{f}(n)$
$\displaystyle:=-\frac{i}{2}\left[(\bar{\psi}\gamma^{5}\psi)_{f,n-1}+(\bar{\psi}\gamma^{5}\psi)_{f,n}\right],$
$\displaystyle=\frac{i}{4a}(-1)^{n}(\chi_{f,n-1}^{\dagger}\chi_{f,n}-\chi_{f,n}^{\dagger}\chi_{f,n-1}-\chi_{f,n}^{\dagger}\chi_{f,n+1}+\chi_{f,n+1}^{\dagger}\chi_{f,n}),$
(49)

for $n=1,2,\cdots,N-2$. Note that both of $S_{f}(n)$ and $PS_{f}(n)$ are
composed of the staggered fermions at the three sites $n$ and $n\pm 1$.

###  3.4 Global observables

We will define the quantum number ($J_{z},\bm{J}^{2},C$, and $P$) and momentum
operators, which will be useful to distinguish the eigenstates of the
Hamiltonian. These operators can be described by some global observables,
which act on the whole lattice.

First of all, let us focus on the isospin operator (4). We define the lattice
version in terms of the staggered fermion. The isospin $J_{z}$ operator counts
the number of particles of each flavor with the factor $\pm 1/2$ on even sites
and the number of anti-particles with the opposite sign on odd sites. Thus, it
can be realized by

$J_{z}=\frac{1}{2}\sum_{n=0}^{N-1}\left(\chi_{1,n}^{\dagger}\chi_{1,n}-\chi_{2,n}^{\dagger}\chi_{2,n}\right).$
(50)

It is convenient to define the isospin $J_{\pm}$ operators by

$J_{\pm}=J_{x}\pm iJ_{y}.$ (51)

Based on the role of fermion operators (26) and (27), $J_{+}$ operator is
given by

$J_{+}=\sum_{n=0}^{N-1}\chi_{1,n}^{\dagger}\chi_{2,n},$ (52)

which transforms $f=2$ particle to $f=1$ particle on even sites and $f=1$
anti-particle to $f=2$ anti-particle on odd sites. Similarly, $J_{-}$ operator
is given by

$J_{-}=\sum_{n=0}^{N-1}\chi_{2,n}^{\dagger}\chi_{1,n},$ (53)

which transforms $f=1$ particle to $f=2$ particle on even sites and $f=2$
anti-particle to $f=1$ anti-particle on odd sites. Then the Casimir operator
$\bm{J}^{2}$ can also be defined as the combination of the operators above by

$\bm{J}^{2}=\frac{1}{2}(J_{+}J_{-}+J_{+}J_{-})+J_{z}^{2}.$ (54)

Second, we will consider the charge conjugation and parity operators. For this
purpose, let us discuss the description of the particle and anti-particle as a
spin state. Applying the Jordan-Winger transformation, the spin representation
of the charge density operator (28) is given by

$\rho_{f,n}=\frac{\sigma_{f,n}^{z}+1}{2}+\frac{(-1)^{n}-1}{2}=\begin{cases}(\sigma_{f,n}^{z}+1)/2&n:\textrm{even},\\\
(\sigma_{f,n}^{z}-1)/2&n:\textrm{odd}.\end{cases}$ (55)

This operator counts the number of particles with $+1$ on even sites and the
number of anti-particles with $-1$ on odd sites. We can confirm that the
particle is described by the spin-up state $\ket{\uparrow}$ on even sites by
taking the expectation value

$\bra{\uparrow}\rho_{f,n}\ket{\uparrow}_{f,n}=\begin{cases}1&n:\textrm{even},\\\
0&n:\textrm{odd}.\end{cases}$ (56)

Similarly, we find that the anti-particle is described by the spin-down state
$\ket{\downarrow}$ on odd sites as

$\bra{\downarrow}\rho_{f,n}\ket{\downarrow}_{f,n}=\begin{cases}0&n:\textrm{even},\\\
-1&n:\textrm{odd}.\end{cases}$ (57)

Based on this fact, charge conjugation, namely the exchange of particles and
anti-particles can be performed by the exchange of even sites and odd sites.
In addition, the spin-up state should be replaced by the spin-down state, and
vice versa. These operations can be realized by the 1-site translation of the
lattice and the multiplication of $\sigma^{x}$ operators. Thus, we define the
charge conjugation operator by Banuls:2013jaa

$C:=\prod_{f=1}^{N_{f}}\left(\prod_{n=0}^{N-1}\sigma_{f,n}^{x}\right)\left(\prod_{n=0}^{N-2}(\mathrm{SWAP})_{f;N-2-n,N-1-n}\right),$
(58)

where the swap operator is given by

$(\mathrm{SWAP})_{f;j,k}=\frac{1}{2}\left(\bm{1}_{f,j}\bm{1}_{f,k}+\sum_{a}\sigma_{f,j}^{a}\sigma_{f,k}^{a}\right),$
(59)

using the Pauli matrices. As the name suggests, the swap operator exchanges
the state $\ket{s}_{f,j}$ and $\ket{s^{\prime}}_{f,k}$, namely

$(\mathrm{SWAP})_{f;j,k}\ket{s}_{f,j}\otimes\ket{s^{\prime}}_{f,k}=\ket{s^{\prime}}_{f,j}\otimes\ket{s}_{f,k}.$
(60)

The product of the swap operators in (58) realizes the 1-site translation. The
charge conjugation defined in this way satisfies $C^{\dagger}C=1$, but
$C^{2}\neq 1$. When we take the periodic boundary condition, $C^{2}=1$ is
achieved in the continuum limit, but this is not the case for the open
boundary condition. Moreover, the Hamiltonian does not commute with $C$ due to
the presence of the boundaries, and we will actually see the eigenstates of
the Hamiltonian give $|\Braket{C}|<1$. Therefore, $C$ does not give a good
quantum number when we take the staggered-fermion regularization with the open
boundary condition. In this study, following the observation of Ref.
Banuls:2013jaa , we assume that the sign of $\mathrm{Re}\Braket{C}$ remembers
the original sign of $C$ for each eigenstate. We discuss this prescription in
detail in Appendix B.

Next, we define the parity operator. The parity transformation
$x\rightarrow-x$ can be achieved by flipping the order of the lattice sites.
The site $n\in\\{0,1,\cdots,N-1\\}$ is mapped to the site $n^{\prime}=N-1-n$.
However, this operation also exchanges particles and anti-particles since the
roles of even sites and odd sites are exchanged when $N$ is even. Thus, an
additional operation of 1-site translation is necessary to fix it. We define
the parity operator by

$\displaystyle P:=\prod_{f=1}^{N_{f}}$
$\displaystyle\left(\prod_{j=0}^{N/2-1}\sigma_{f,2j+1}^{z}\right)$
$\displaystyle\times$
$\displaystyle\left(\prod_{n=0}^{N-2}(\mathrm{SWAP})_{f;N-2-n,N-1-n}\right)\left(\prod_{n=0}^{N/2-1}(\mathrm{SWAP})_{f;n,N-1-n}\right),$
(61)

where the products of the swap operators perform the reversal $n\rightarrow
n^{\prime}$ and the 1-site translation.333 If we implement the reversal
$n\rightarrow n^{\prime}$ in this manner, the bond dimension of MPO grows
exponentially with $N$. Thus, in practice, we apply the reversal by
transposing all the matrices in MPS. The additional factors of $\sigma^{z}$
come from the shift of the staggered phase, which corresponds to $\gamma^{0}$
in the parity transformation of the Dirac fermion
$\psi(x)\rightarrow\gamma^{0}\psi(-x)$. As we mentioned for the $C$ operator,
the $P$ operator in the open boundary condition does not commute with the
Hamiltonian as it contains the $1$-unit lattice translation. Therefore, we
take the same prescription to determine the parity quantum number for each
state as in the case of $C$.

Finally, the other important quantity is a total momentum operator, which can
be used to identify the momentum excitation Banuls:2013jaa . We start with the
continuum description of the gauge invariant operator

$K=\sum_{f=1}^{N_{f}}\int
dx\,\psi_{f}^{\dagger}(i\partial_{x}-A_{1})\psi_{f},$ (62)

which commutes with the continuum Hamiltonian (20) under the periodic boundary
condition using the Gauss-law constraint (21). In our case with the open
boundary condition, it does not commute with the Hamiltonian since the
translational symmetry is explicitly broken. Thus, the expectation value
$\Braket{K}$ is no longer the quantum number in the strictest sense. However,
we will see that the operator is still useful as an approximate one to
investigate the mass spectrum of the model. Let us consider its lattice
version. Here we set $A_{1}(x)=0$ since we fix the gauge $U_{n}=1$ in our
setup. The combination $\psi_{f}^{\dagger}\partial_{x}\psi_{f}$ of the Dirac
fermion corresponds to

$\psi_{f}^{\dagger}\partial_{x}\psi_{f}=(\psi_{u}^{\dagger}\partial_{x}\psi_{u}+\psi_{d}^{\dagger}\partial_{x}\psi_{d})_{f}=\frac{1}{2a}\chi_{f,n}^{\dagger}(\chi_{f,n+2}-\chi_{f,n}),$
(63)

in terms of the staggered fermion. There is another possible combination

$-(\partial_{x}\psi_{f}^{\dagger})\psi_{f}=-\frac{1}{2a}(\chi_{f,n+2}^{\dagger}-\chi_{f,n}^{\dagger})\chi_{f,n},$
(64)

given by the integral by parts ignoring boundary term. Then we define the
total momentum on the lattice as a Hermitian operator by taking symmetric
combination

$\displaystyle K$
$\displaystyle:=\frac{i}{2}\sum_{f=1}^{N_{f}}\sum_{n=0}^{N-3}\frac{1}{2a}\left[\chi_{f,n}^{\dagger}(\chi_{f,n+2}-\chi_{f,n})-(\chi_{f,n+2}^{\dagger}-\chi_{f,n}^{\dagger})\chi_{f,n}\right],$
$\displaystyle=\frac{i}{4a}\sum_{f=1}^{N_{f}}\sum_{n=1}^{N-2}(\chi_{f,n-1}^{\dagger}\chi_{f,n+1}-\chi_{f,n+1}^{\dagger}\chi_{f,n-1}).$
(65)

This operator does not commute with the term $H_{w}$ (34) and $H_{J}$ (33) of
the lattice Hamiltonian due to the open boundary. We also note that the latter
$[K,H_{J}]$ has an $O(a)$ violation effect even in the periodic boundary
condition.

##  4 Calculation method and the simulation setup

We employ the density matrix renormalization group (DMRG) White:1992zz ;
White:1993zza ; Schollw_ck_2005 ; Schollw_ck_2011  to study the spin
Hamiltonian of the 2-flavor Schwinger model after the Jordan-Winger
transformation, whose explicit form is given by (95). The DMRG is known as an
efficient method to study (1+1)d gapped spin systems and has been developed
mainly in the field of condensed matter physics. We utilized the C++ library
of ITensor itensor  to perform the tensor network calculation of this work.
Let us briefly explain the basic idea of DMRG to obtain the ground state and
excited states, and then we explain the details of parameter settings.

###  4.1 Quick review of DMRG

In the spin systems, any wave function $\ket{\Psi}$ can be expressed as the
form of the matrix product states (MPS),

$\ket{\Psi}=\sum_{i_{1},\ldots,i_{N}=1}^{2}\tr[A_{1}^{(i_{1})}\cdots
A_{N}^{(i_{N})}]\ket{i_{1}\ldots i_{N}},$ (66)

by repeating the singular-value decomposition (SVD). Here, $i_{n}=1,2$ denotes
the spin degrees of freedom at the $n$-th site, $A_{n}^{(i_{n})}$ denotes a
$D\times D$ matrix, and this size $D$ is called the bond dimension. The upper
bound for the entanglement entropy of $\ket{\Psi}$ is given by $\ln D$.
Therefore, the MPS gives a useful tool to study the many-body states with low
entanglement entropies, such as the ground state of the $(1+1)$d gapped
systems Stoudenmire_2012 ; Wall_2012 . For the 2-flavor Schwinger model, we
arrange the site index $n$ and the flavor index $f$ on the 1d lattice with the
single index to apply DMRG. The ordering of the indices is chosen so that the
behavior of entanglement entropy is reproduced appropriately with a reasonable
bond dimension. This point is discussed in Appendix C.

The DMRG is a variational algorithm based on the MPS. In each step of the
algorithm, the matrices are updated to decrease the energy
$E=\Braket{\Psi}{H}{\Psi}$ as a cost function. In addition, we perform the
low-rank approximation and thus the smaller singular values are discarded,
which amount to an error $\Delta$. We determine the bound dimension by setting
the maximal bond dimension and also the cutoff parameter $\varepsilon$ on the
error so that $\Delta\leq\varepsilon$. Smaller $\varepsilon$ gives a better
approximation, but it also requires a larger bond dimension and increases the
computational costs. We can also effectively calculate the expectation values
or correlation functions of local operators by rewriting those operators in
the form of matrix product operators (MPOs) and then taking contractions with
the ground state $\ket{\Psi}$.

We can use DMRG to obtain the low-energy excited states in a recursive way.
Assume that we already find the energy eigenstates $\ket{\Psi_{k^{\prime}}}$
with $k^{\prime}=0,1,\ldots,k-1$ from below. Then, we apply the same technique
to find the $k$-th state $\ket{\Psi_{k}}$ by changing the Hamiltonian for the
cost function as

$H_{k}=H+W\sum_{k^{\prime}=0}^{k-1}\ket{\Psi_{k^{\prime}}}\bra{\Psi_{k^{\prime}}},$
(67)

where $W>0$ is a weight to impose the orthogonality. We can generate the
excited states from the ground state to any level step by step.

###  4.2 Simulation setup

Let us explain our parameter setup when using the ITensor itensor . The gauge
coupling $g$ has mass dimension $1$ in $1+1$d QED, and thus we can measure the
energy scale in the unit of $g$ by setting $g=1$. In this work, we always set
$g=1$ and the fermion mass $m=0.1$, so the photon mass is
$\mu=\sqrt{\frac{2}{\pi}}\simeq 0.8$. The lattice fermion mass (31) becomes
$m_{\mathrm{lat}}=0.1-\frac{a}{4}$. The theta angle is normally set to
$\theta=0$, except when measuring the one-point function of the pion at
$\theta=2\pi$.

For the correlation-function scheme and the one-point-function scheme, we use
the lattice size of $N=160$. The lattice spacing is set to $a\approx 0.25$ so
that the physical size is $L=a(N-1)=39.8$. The number of DMRG steps called the
sweeps, is set to $N_{\mathrm{sweep}}=20$. We generate the ground state for
four different values of the cutoff parameter: $\varepsilon=10^{-10}$,
$10^{-12}$, $10^{-14}$, and $10^{-16}$. To characterize the bond dimension of
the MPS, we focus on the largest number of nonzero singular values, which we
call the effective bond dimension, denoted as $D_{\mathrm{eff}}$. In our
computations, we set the maximal bond dimension large enough so that
$D_{\mathrm{eff}}$ is solely controlled by the cutoff $\varepsilon$ for the
above physical setup. We observe $D_{\mathrm{eff}}$ to be about 400, 800,
1600, and 2800 for the respective values of $\varepsilon$ above.

For the dispersion-relation scheme, we generate many excited states up to
$k=23$, which require a lot of computational cost. Therefore, we choose a
smaller lattice size of $L=19.8$ with $N=100$ and $a=0.2$. The excited states
are generated with a cutoff of $\varepsilon=10^{-10}$ and a weight parameter
of $W=10$. To achieve better convergence of higher states, we increase the
number of sweeps to $N_{\mathrm{sweep}}=50$. The bond dimension is about 500
for the ground state while it is at most 2300 for the excited states.

As an initial state of the DMRG, we choose the Néel state, which is a direct
product of the spin-down states on even sites and the spin-up state on odd
sites,

$\Ket{\mathrm{N\acute{e}el}}=\bigotimes_{f=1}^{N_{f}}\ket{\downarrow}_{f,0}\ket{\uparrow}_{f,1}\cdots\ket{\downarrow}_{f,N-2}\ket{\uparrow}_{f,N-1}.$
(68)

Based on (56) and (57), the Néel states is regarded as a zero-particle state.
We also impose the charge conservation condition during the DMRG, so that the
MPS satisfies the condition $Q=0$, where

$Q=\sum_{f=1}^{N_{f}}\sum_{n=0}^{N-1}\rho_{f,n}=\frac{1}{2}\sum_{f=1}^{N_{f}}\sum_{n=0}^{N-1}\sigma_{f,n}^{z}.$
(69)

We note that the Gauss law with the usual open boundary on both sides requires
$Q=0$ on the physical states.

##  5 Simualtion results

In this section, we explain our numerical results for the meson spectrum of
the $2$-flavor Schwinger model at $\theta=0$. We apply three distinct methods
in our computations of the meson spectrum:

  * •

the correlation-function scheme,

  * •

the one-point-function scheme, and

  * •

the dispersion-relation scheme.

Each method has its own pros and cons, and we are going to discuss them. We
will see that all these schemes give consistent results.

###  5.1 Correlation-function scheme

In the relativistic quantum field theories, the Hilbert space only plays the
secondary role, and we are supposed to reconstruct all the physical
information from the correlation functions. In the conventional Euclidean
lattice gauge theory, people usually follow this dogma, and the mass spectrum
is obtained from the correlation function in the imaginary time direction. We
can take the same approach also in the Hamiltonian formalism by the
measurement of the spatial correlation function with the distance $r=|x-y|$.

First, let us work on pions, and we consider the equal-time spatial
correlation function,

$C_{\pi}(r)=\Braket{\pi(x)\pi(y)},$ (70)

where $\pi(x)$ denotes the operator defined by (43) with $x=na$. In order to
evade the boundary effect as much as possible, we compute $C_{\pi}(r)$ by
changing $x$ and $y$ symmetrically as $x=(L-r)/2$ and $y=(L+r)/2$, and the
range of $r$ is restricted to $0\leq r\leq L/2$. The results are shown in the
left panel of Fig. 1 in the logarithmic scale, and the pion mass can be
extracted from the exponential decay of $C_{\pi}(r)$. Here, the data with
different colors represent the different values of the cutoff parameter
$\varepsilon$.

Figure 1:  (Left) The correlation function of pion
$\ln|\Braket{\pi(x)\pi(y)}|$ is plotted against the distance $r=|x-y|$ for
various values of $\varepsilon$. The number of lattice sites is $N=160$ and
the lattice spacing $a$ is determined so that $L=a(N-1)=39.8$. (Right) The
effective mass of the pion $M_{\pi,\mathrm{eff}}(r)$ (3-point average)
calculated from the correlation function in the left panel is plotted against
$r$.

It is convenient to use the so-called effective mass defined by

$\tilde{M}_{\pi,\mathrm{eff}}(r)=-\frac{1}{2a}\log\frac{C_{\pi}(r+2a)}{C_{\pi}(r)},$
(71)

where $2a$ comes from the step size of changing $r$. We further take the
3-point average of the effective mass

$M_{\pi,\mathrm{eff}}(r)=\frac{1}{4}\tilde{M}_{\pi,\mathrm{eff}}(r-2a)+\frac{1}{2}\tilde{M}_{\pi,\mathrm{eff}}(r)+\frac{1}{4}\tilde{M}_{\pi,\mathrm{eff}}(r+2a)$
(72)

to suppress a remaining oscillation caused by the staggered fermion. The
result of $M_{\pi,\mathrm{eff}}(r)$ is shown in the right panel of Fig. 1.

One might be tempted to think that the pion mass corresponds to the plateau
value of $M_{\pi,\mathrm{eff}}(r)$, and the result with the DMRG cutoff
$\varepsilon=10^{-10}$ actually seems to become constant for $r\gtrsim 10$
almost exactly. However, this is the fake plateau due to the low-rank
approximation. The point is that the leading asymptotic behavior of the
spatial correlator is not purely the exponential decay, and it would take the
Yukawa-type form asymptotically as $r\to\infty$,

$C_{\pi}(r)\sim\frac{1}{r^{\alpha}}\exp(-M_{\pi}r).$ (73)

We actually have $\alpha=1/2$ for the $(1+1)$d free massive boson, and we
shall discuss the detailed analysis in Appendix D in the case of the
$1$-flavor Schwinger model. As a result, the effective mass for the Yukawa-
type correlation function is given by

$M_{\pi,\mathrm{eff}}(r)=-\frac{d}{dr}\log
C_{\pi}(r)\sim\frac{\alpha}{r}+M_{\pi},$ (74)

and there must be an additional $O(1/r)$ contribution on top of the actual
mass $M_{\pi}$.

Figure 2:  The effective mass of the pion $M_{\pi,\mathrm{eff}}(r)$ is plotted
against $1/r$. The data points for $\varepsilon=10^{-16}$ are fitted by
$\alpha/r+M$ inside the region $0.075\leq 1/r\leq 0.15$. The fitting result is
depicted by the shaded band with systematic error.

Motivated by this fact, we plot $M_{\pi,\mathrm{eff}}(r)$ against $1/r$ in
Fig. 2. We can see that the behavior of the effective mass strongly depends on
the cutoff $\varepsilon$ especially when $r$ is large. When $\varepsilon$ is
not small enough, we observe the saturation of $M_{\pi,\mathrm{eff}}(r)$, and
then the $\alpha/r$ term seems to be absent. We note that the low-rank
approximation of the DMRG is similar to the approximation of the transfer
matrix by a finite matrix. Thus, $C_{\pi}(r)$ in the DMRG is approximated by
the sum of purely exponential functions, and we need sufficiently large bond
dimensions to reproduce the non-exponential corrections, such as
$1/r^{\alpha}$.

In fact, we can observe in Fig. 2 that the development of the $1/r$-behavior
in $M_{\pi,\mathrm{eff}}(r)$ for large $r$ by making $\varepsilon$
sufficiently small, i.e. the bond dimension sufficiently large. We estimate
the mass $M_{\pi}$ by the linear extrapolation $1/r\rightarrow 0$ of the
result for the largest bond dimension with $\varepsilon=10^{-16}$, which is
performed by fitting the data points with $\alpha/r+M_{\pi}$. To evaluate the
systematic errors from the uncertainty of the fitting range, we try a lot of
fittings by changing the fitting range inside the region $0.075\leq 1/r\leq
0.15$, and we obtain the probability distribution of the fitting results. The
best-fitting result and its error are estimated from the position and the
width of the peak, respectively. Thus, we obtained

$M_{\pi}=0.431(1),$ (75)

with $\alpha=0.477(9)$, and the fitting lines are drawn as the purple shadow
in Fig. 2.

Figure 3:  The effective mass of sigma meson $M_{\sigma,\mathrm{eff}}(r)$
(left) and eta meson $M_{\eta,\mathrm{eff}}(r)$.

Next, we perform similar analyses for sigma meson (45) and eta meson (44).
Since these are isospin singlets, their one-point functions are not zero, and
we subtract the disconnected parts from the correlation functions,

$\displaystyle C_{\sigma,\mathrm{eff}}(r)$
$\displaystyle=\langle\sigma(x)\sigma(y)\rangle-\langle\sigma(x)\rangle\langle\sigma(y)\rangle,$
(76) $\displaystyle C_{\eta,\mathrm{eff}}(r)$
$\displaystyle=\langle\eta(x)\eta(y)\rangle-\langle\eta(x)\rangle\langle\eta(y)\rangle,$
(77)

with $x=(L-r)/2$, $y=(L+r)/2$. We then compute the 3-point averages of the
effective mass, $M_{\sigma,\mathrm{eff}}(r)$ and $M_{\eta,\mathrm{eff}}(r)$,
and they are shown in Fig. 3. The difference in the asymptotic behavior is
observed by changing $\varepsilon$ also in these cases. We plot the effective
masses of the sigma and eta mesons against $1/r$ in Fig. 4 to see the
asymptotic behavior. They approach $\propto 1/r$ for smaller $\varepsilon$ as
expected. We fit the data for $\varepsilon=10^{-16}$ by $\alpha/r+M$ inside
the region $0.075\leq 1/r\leq 0.15$ and estimate the systematic error. Then we
obtained

$M_{\sigma}=0.722(6),$ (78)

with $\alpha=0.83(5)$ for sigma meson, and

$M_{\eta}=0.899(2),$ (79)

with $\alpha=0.51(2)$ for eta meson. It is notable that $\alpha_{\sigma}\sim
0.8$ has a relatively large deviation from the free boson result,
$\alpha=1/2$, which may suggest that the sigma meson has a nontrivial
dispersion relation even for small momentum.

Figure 4:  The effective mass of sigma meson $M_{\sigma,\mathrm{eff}}(r)$
(left) and eta meson $M_{\eta,\mathrm{eff}}(r)$ (right) are plotted against
$1/r$. The data points for $\varepsilon=10^{-16}$ are fitted by $\alpha/r+M$
inside the range $0.075\leq 1/r\leq 0.15$. The fitting results are depicted by
the shaded bands with systematic errors.

Finally, we summarize the masses of the three mesons measured by the
correlation functions:

$\begin{array}[]{c|c|c|c}&\text{pion}&\text{sigma}&\text{eta}\\\
\hline\cr\text{mass}/g&\,\,0.431(1)&\,\,0.722(6)&\,\,0.899(2)\end{array}$ (80)

The numerical results are qualitatively consistent with the analytic result by
bosonization $M_{\pi}<M_{\sigma}<M_{\eta}$. We also find

$M_{\sigma}/M_{\pi}=1.68(2),$ (81)

which is close to the prediction by the sine-Gordon model Eq. (19).

###  5.2 One-point-function scheme

We consider an alternative way to obtain the mass spectrum without using the
two-point correlation functions. Let us recall that we are taking the open
boundary condition, and we can use those boundaries as the source for
excitations from the thermodynamic ground state. The boundary effect decays
exponentially for the gapped systems, and thus the one-point function of a
local operator $\mathcal{O}(x)$ should behave as
$\langle\mathcal{O}\rangle+C^{\prime}e^{-M_{\mathcal{O}}x}$ as the function of
the distance $x=an$ from the boundary. Here, $\langle\mathcal{O}\rangle$ gives
the vacuum expectation value in the thermodynamic limit, and $M_{\mathcal{O}}$
in the exponent gives the lightest particle mass with the same quantum number
of $\mathcal{O}(x)$. In the context of condensed matter physics, it is known
that the correlation function can be obtained from the Fridel oscillation,
which is induced by a boundary effect or a local external field
PhysRevB.54.13495 ; SHIBATA19971024 .

We note that the $x$-dependence in this method takes the purely exponential
form $e^{-Mx}$ as the leading behavior for $x\to\infty$. This can be easily
understood by considering the path integral and the $\pi/2$ rotation of
Euclidean spacetime. Then, the boundary condition sits at the constant
imaginary time and defines the state $\ket{\mathrm{Bdry}}$ with zero momentum.
Thus, the leading contribution to the imaginary-time correlation function
$\bra{\mathrm{Vac}}Oe^{-H|x|}\ket{\mathrm{Bdry}}$ should come from the
lightest particle with the zero-momentum projection, giving $e^{-Mx}$. This
feature has nice compatibility with the low-rank approximation of DMRG.

####  5.2.1 The one-point functions of $\sigma$ and $\eta$ at $\theta=0$

At $\theta=0$, the boundary condition turns out to be completely invariant
under the isospin rotation, and thus the boundary state $\ket{\mathrm{Bdry}}$
does not produce one pion states. Therefore, let us here focus on the iso-
singlet particles, $\sigma$ and $\eta$, and we will come back to pions later.

First, we discuss the eta meson as it turns out to be the simplest one. Since
the $G$ parity is not spontaneously broken, we must have
$\langle\eta\rangle=0$ in the thermodynamic limit. However, the staggered
fermion realizes the $G$ parity (or charge conjugation) as the one-unit
lattice translation, and thus the open boundary condition violates the $G$
parity. Therefore, the boundary state can be a source of the eta meson, and we
evaluate the one-point function $\Braket{\eta(x)}$ of the eta meson operator
(44) in the range $0<x\leq L/2$. The result is shown in Fig. 5. The cutoff
parameter is changed from $\varepsilon=10^{-10}$ to $10^{-16}$. The one-point
function decays exponentially with $x$ as expected. Thus, we fit the data
points of $\ln|\Braket{\eta(x)}|$ by $-M_{\eta}x+C$ in the fitting range
$7\leq x\leq 13$, and the result is

$M_{\eta}=0.9014(1),$ (82)

with $C=-1.096(1)$ for the smallest cutoff $\varepsilon=10^{-16}$. The
corresponding fitting curve is shown in Fig. 5 with the purple line. In this
case, we also find that the results for the other values of $\varepsilon$ are
consistent within the fitting error. Thus, the cutoff dependence does not
appear unlike the case of the correlation functions, and we suspect that this
is because MPS can efficiently express purely exponential decay.

Figure 5:  The one-point function $\ln|\Braket{\eta(x)}|$ of the eta meson is
plotted against $x=an$ with $n=1,\cdots,N/2-1$ for various values of
$\varepsilon$. The number of lattice sites is $N=160$ and the lattice spacing
$a$ is determined so that $L=a(N-1)=39.8$. The result of fitting by
$-M_{\eta}x+C$ for $\varepsilon=10^{-16}$ is also plotted by the solid line
inside the range and by the broken line outside.

Next, we evaluate the one-point function $\Braket{\sigma(x)}$ of the sigma
meson (45) for $0<x\leq L/2$. We note that $\sigma$ has the same quantum
number with the vacuum, and then $\Braket{\sigma(x)}$ is nonzero also in the
bulk. It behaves as $e^{-Mx+C}+A$ with a constant shift of $A$, so we subtract
the value $\Braket{\sigma(L/2)}$ at the center $x=L/2$ of the lattice from
$\Braket{\sigma(x)}$ to remove the constant. The result is shown in Fig. 6,
which indicates the exponential decay as expected. We fit the data points of
$\ln|\Braket{\sigma(x)-\sigma(L/2)}|$ by $-M_{\sigma}x+C$ in the range $7\leq
x\leq 13$, and the best-fit parameter is

$M_{\sigma}=0.761(2),$ (83)

with $C=-2.71(2)$, which are independent of the value of $\varepsilon$. The
result of fitting for $\varepsilon=10^{-16}$ is shown in Fig. 6 with the
purple line.

Figure 6:  The one-point function $\ln|\Braket{\sigma(x)-\sigma(L/2)}|$ of the
sigma meson is plotted against $x=an$ with $n=1,\cdots,N/2-1$ for various
values of $\varepsilon$. The value at $x=L/2$ is subtracted from
$\Braket{\sigma(x)}$ to eliminate the constant shift in the bulk. The number
of lattice sites is $N=160$ and the lattice spacing $a$ is determined so that
$L=a(N-1)=39.8$. The result of fitting by $-M_{\sigma}x+C$ is also plotted by
the solid line inside the range and by the broken line outside.

####  5.2.2 The one-point functions of $\pi_{3}$ at $\theta=2\pi$

Let us now come back to the issue of pions. As we have argued, the boundary
state at $\theta=0$ is neutral under the isospin rotation, and thus it does
not produce one-pion states and we have $\langle\pi(x)\rangle=0$ for all $x$.
Therefore, we need to somehow create the boundary state that transforms
nontrivially under the isospin rotation to study pions with the one-point-
function scheme.

In this study, we decided to use one of the ground states at $\theta=2\pi$ for
this purpose. Since the Hamiltonians at $\theta=0$ and $\theta=2\pi$ are
unitary equivalent under the periodic boundary condition, the bulk properties
are the exactly same between $\theta=0,2\pi$. As we have discussed in Section
2.2, the ground state at $\theta=2\pi$ is a nontrivial SPT state protected by
the isospin $\mathrm{SU}(2)_{V}/\mathbb{Z}_{2}$ symmetry, and thus the
boundary states with the open boundary condition have the isospin $1/2$. This
boundary charge can be a source of the pions so that $\Braket{\pi(x)}$ becomes
nonzero.

About the computational cost, it turns out that the bond dimensions for the
MPS are mostly the same at $\theta=0$ and $\theta=2\pi$ when the system size
is large enough. Therefore, we can obtain the ground state at $\theta=2\pi$ as
easily as that of $\theta=0$. We, however, observe that the bond dimension at
$\theta=2\pi$ increases significantly if the system size is not large enough,
and we suspect its reason is as follows. At $\theta=2\pi$, there is $4$-fold
degeneracy due to the boundary degrees of freedom, but they split into the
singlet and the triplet states with the energy splitting $\sim e^{-M_{\pi}L}$.
That is, the true ground state at finite $L$ has an extra Bell pair between
the boundary isospin $1/2$ states, which adds $\ln 2$ to the entanglement
entropy. When we cut the system at $x=L/2$, this extra $\ln 2$ should be
accumulated by the large numbers of small singular values, and thus the bond
dimension becomes quite huge just to create the Bell pair between the
boundaries. If $L$ is large enough, the energy gain by creating the Bell pair
becomes negligible, and the DMRG would produce one of the ground states with
disentangled boundary states practically. Thus, the computational cost becomes
almost the same as that for the trivial state at $\theta=0$.

Let us now evaluate the one-point function $\Braket{\pi(x)}$ of the pion (43)
for $0<x\leq L/2$ using the ground state at $\theta=2\pi$. The result is shown
in Fig. 7. We again find the exponential decay, and thus fit the data points
of $\ln|\Braket{\pi(x)}|$ by $-M_{\pi}x+C$ in the range $7\leq x\leq 13$. The
result is

$M_{\pi}=0.4175(9),$ (84)

with $C=0.203(9)$, which do not depend on the cutoff $\varepsilon$ as before.
The fitting result for $\varepsilon=10^{-16}$ is shown in Fig. 7 with the
purple line.

Figure 7:  The one-point function $\ln|\Braket{\pi(x)}|$ of the pion is
plotted against $x=an$ with $n=1,\cdots,N/2-1$ for various values of
$\varepsilon$. The number of lattice sites is $N=160$ and the lattice spacing
$a$ is determined so that $L=a(N-1)=39.8$. We set $\theta=2\pi$ in order to
induce the boundary charges, which make $\Braket{\pi(x)}$ nonzero. The result
of fitting by $-M_{\pi}x+C$ is also plotted by the solid line inside the range
and by the broken line outside.

Let us summarize the effective masses using the one-point-function scheme:

$\begin{array}[]{c|c|c|c}&\text{pion}&\text{sigma}&\text{eta}\\\
\hline\cr\text{mass}/g&\,\,0.4175(9)&\,\,0.761(2)&\,\,0.9014(1)\end{array}$
(85)

The order of three meson masses is consistent with the analytic prediction. We
also find

$M_{\sigma}/M_{\pi}\simeq 1.821(6),$ (86)

which is still close to the WKB prediction, $\sqrt{3}$, with a 5% deviation.
The significant feature of the method is that the results do not depend on the
cutoff parameter $\varepsilon$ as long as it is sufficiently small. Therefore,
the systematic error from the cutoff is expected to be small enough. We do not
need to increase the bond dimension so much, unlike the method by the
correlation function.

###  5.3 Disparsion-relation shceme

So far we have studied the mass spectrum by using the local observables of the
ground state, and these methods are applicable both in the path integral and
the Hamiltonian formalisms. As the third method for computing the mass
spectrum, let us take a different approach that is specific to the Hamiltonian
formalism: We compute the excited states as explained in Section 4.1, and then
determine the mass spectrum from the dispersion relation.

The low-lying excited states correspond to one-particle excitations. For
example, the zero modes of the lightest meson, namely the pion, is expected to
be obtained as the first excited state. We can also obtain the states with
nonzero momentum $K$, and we can fit the data with the dispersion relation
$\Delta E\simeq\sqrt{K^{2}+M_{\pi}^{2}}$ to obtain the pion mass. As we go to
the higher excited states, we will encounter one-particle states of the sigma
and eta mesons. They can be distinguished by measuring quantum numbers, such
as the isospin and $G$-parity. Thus, we can compute the mass spectrum from the
dispersion relation by generating the excited states.

We note that our computation is done in the finite open interval, and thus the
momentum is not a good quantum number. Also, there may exist a nontrivial
contribution to the excitation energy from the boundaries. We are neglecting
those subtleties in this work, but, surprisingly, it turns out that the
numerical results are almost consistent with those with the previous two
methods.

We generated the MPS up to the 23rd excited state at $\theta=0$ with the small
physical volume $L=19.8$. The energy gap $\Delta E_{k}=E_{k}-E_{0}$ of the
$k$-th excited state is shown in the left panel of Fig. 8. We also measured
the square of total momentum $K^{2}$ defined by (65). We note that its ground-
state expectation value $\Braket{K^{2}}_{0}\simeq 0.46$ is nonzero because of
the boundary effect and maybe also due to lattice artifacts, and thus we
subtract $\Braket{K^{2}}_{0}$ from $\Braket{K^{2}}_{k}$ of the excited states.
The result is plotted in the right panel of Fig. 8. From these results, we
find many triply degenerated states, which are candidates for the states of
the pion. There are a few singlet states as well, which are candidates for the
eta and sigma mesons.

Figure 8:  (Left) The energy gap $\Delta E_{k}=E_{k}-E_{0}$ is plotted against
the level of the excited state $k$. (Right) The square of total momentum
$\Delta K_{k}^{2}=\Braket{K^{2}}_{k}-\Braket{K^{2}}_{0}$ is plotted against
$k$ after subtracting the result for the ground state.

To identify the states, we measure the expectation values of the isospin
operators, $\bm{J}^{2}$ and $J_{z}$, the parity $P$ and the $G$-parity
$G=Ce^{i\pi J_{y}}$ in Section 3.4. We note that the DMRG does not produce the
states in a diagonal basis for these quantities. We diagonalize the $3\times
3$ matrix $\Braket{\psi_{k_{1}}}{J_{z}}{\psi_{k_{2}}}$ in each triplet to
compute expectation values in the $J_{z}$ basis.444In computing the
expectation value of the $G$-parity, we find it easier to do it in the $J_{y}$
basis instead of the $J_{z}$ basis because $G=Ce^{i\pi J_{y}}$, and we thus
performed it in the $J_{y}$ basis.,555It is possible that triplets and
singlets are also mixed in the DMRG if their energies are close. In fact, the
states for $k=19,\cdots,23$ are mostly degenerated. We separated one triplet
and two singlets out of them by diagonalization of
$\Braket{\psi_{k_{1}}}{\bm{J}^{2}}{\psi_{k_{2}}}$ and of
$\Braket{\psi_{k_{1}}}{C}{\psi_{k_{2}}}$.

$k$ | $\bm{J}^{2}$ | $J_{z}$ | $G$ | $P$
---|---|---|---|---
1 | 2.00000004 | 0.99999997 | 0.27872443 | -6.819$\times{10}^{-8}$
2 | 2.00000012 | -0.00000000 | 0.27872416 | -6.819$\times{10}^{-8}$
3 | 2.00000004 | -0.99999996 | 0.27872443 | -6.819$\times{10}^{-8}$
4 | 2.00000007 | 0.99999999 | 0.27736066 | 7.850$\times{10}^{-8}$
5 | 2.00000006 | 0.00000000 | 0.27736104 | 7.850$\times{10}^{-8}$
6 | 2.00000009 | -0.99999998 | 0.27736066 | 7.850$\times{10}^{-8}$
7 | 2.00000010 | 1.00000000 | 0.27536687 | -8.838$\times{10}^{-8}$
8 | 2.00000002 | 0.00000000 | 0.27536702 | -8.837$\times{10}^{-8}$
9 | 2.00000007 | -0.99999998 | 0.27536687 | -8.838$\times{10}^{-8}$
10 | 2.00000007 | 0.99999998 | 0.27356274 | 9.856$\times{10}^{-8}$
11 | 2.00000005 | 0.00000001 | 0.27356277 | 9.856$\times{10}^{-8}$
12 | 2.00000007 | -0.99999999 | 0.27356274 | 9.856$\times{10}^{-8}$
15 | 1.99999942 | 0.99999966 | 0.27173470 | -1.077$\times{10}^{-7}$
16 | 2.00000052 | 0.00000000 | 0.27173482 | -1.077$\times{10}^{-7}$
17 | 2.00000015 | -1.00000003 | 0.27173470 | -1.077$\times{10}^{-7}$
19 | 2.00009067 | 1.00004377 | 0.27717104 | -3.022$\times{10}^{-8}$
20 | 2.00002578 | -0.00000004 | 0.27717020 | -3.023$\times{10}^{-8}$
21 | 2.00003465 | -1.00001622 | 0.27717104 | -3.023$\times{10}^{-8}$
Table 1:  The quantum numbers of the isospin triplet states. The index $k$ comes from the level of each state in the original basis. The rows of the table are separated into each triplet. $k$ | $\bm{J}^{2}$ | $J_{z}$ | $G$ | $P$
---|---|---|---|---
0 | 0.00000003 | -0.00000000 | 0.27984227 | 3.896$\times{10}^{-7}$
13 | 0.00000003 | 0.00000000 | 0.27865844 | 1.273$\times{10}^{-7}$
14 | 0.00000003 | 0.00000000 | 0.27508176 | -2.765$\times{10}^{-8}$
18 | 0.00000028 | 0.00000006 | -0.27390909 | -6.372$\times{10}^{-7}$
22 | 0.00001537 | 0.00000115 | 0.26678987 | 7.990$\times{10}^{-8}$
23 | 0.00003607 | -0.00000482 | -0.27664779 | 5.715$\times{10}^{-7}$
Table 2:  The quantum numbers of the isospin singlet states.

The expectation values of $\bm{J}^{2}$, $J_{z}$, $G$, and $P$ in the $J_{z}$
basis are listed in Tables 1 and 2 for iso-triplets and iso-singlets,
respectively. The index $k$ comes from the level of the state on the original
random basis. We find that $|G|\neq 1$ because of $|C|\neq 1$ by the effect of
the boundary. Hopefully, the sign of $G$ can be assumed to remember the
original quantum number Banuls:2013jaa , and, if it is true, we can still
identify the $G$-parity. This point will be discussed more in details in
Appendix B. We identify the lowest triplet $k=1,2,3$ as the lowest modes of
the pions ($\pi^{+}$, $\pi^{0}$, $\pi^{-}$) since they have the quantum
numbers consistent with the pion, namely $J^{PG}=1^{-+}$ and $J_{z}=0,\pm 1$.
For the iso-singlets shown in Table 2, we find that the $k=13$ state has the
quantum number consistent with the sigma meson, namely $J^{PG}=0^{++}$ and
$J_{z}=0$. The $k=18$ state is consistent with the eta meson with
$J^{PG}=0^{--}$ and $J_{z}=0$. We identify these singlets with the lowest
modes of the sigma and eta mesons.

Figure 9:  The energy gap $\Delta E_{k}$ is plotted against the square of
total momentum $\Delta K_{k}^{2}$. The states with the same isospin and
$G$-parity are plotted by the same symbol. Then each state is identified with
the pion, sigma, or eta meson. We fit the data for each meson by $\Delta
E=\sqrt{b^{2}\Delta K^{2}+M^{2}}$. The results are shown by the broken lines.
The values of $M$ for each meson are also plotted as the endpoints of the
fitting lines.

After identifying the quantum numbers, we plot the energy gap $\Delta
E_{k}=E_{k}-E_{0}$ against the momentum square $\Delta
K_{k}^{2}=\Braket{K^{2}}_{k}-\Braket{K^{2}}_{0}$ to obtain the dispersion
relation as shown in Fig. 9. The states with the same isospin $\bm{J}^{2}$ and
$G$-parity are plotted by the same symbol.666 The triplet for $k=19,20,21$ is
not shown in this plot since it is not of the state of the single pion. We
expect that the triplet comes from the pion scattering state, which was
discussed in Harada:1993va . Then we fit the data points by $\Delta
E=\sqrt{b^{2}\Delta K^{2}+M^{2}}$ with fitting parameters $M$ and $b$. The
fitting result of $M$ can be regarded as the mass of the corresponding meson
as an extrapolation to $\Delta K^{2}\rightarrow 0$. We obtained
$M_{\pi}=0.426(2)$, $b_{\pi}=1.017(4)$ for the pion; and
$M_{\sigma}=0.7456(5)$, $b_{\sigma}=1.087(2)$ for the sigma meson. The fitting
for the eta meson is simply solving an equation since there are only two data
points. The result are $M_{\eta}=0.904$ and $b_{\eta}=0.962$. We summarize the
masses of the mesons determined by the energy gap of the excited states:

$\begin{array}[]{c|c|c|c}&\text{pion}&\text{sigma}&\text{eta}\\\
\hline\cr\text{mass}/g&\,\,0.426(2)&\,\,0.7456(5)&\,\,0.904\end{array}$ (87)

We find the mass ratio

$M_{\sigma}/M_{\pi}\simeq 1.75(1)$ (88)

from this result, which is close to the WKB prediction $\sqrt{3}$.

##  6 Conclusion and Discussion

In this paper, we work on three independent methods to compute the mass
spectrum of lattice gauge theories in the Hamiltonian formalism, which apply
to tensor networks and quantum computation. The methods are tested in the
massive 2-flavor Schwinger model at $\theta=0$, some of which properties are
analogous to the ones of $4$d QCD. The two species of fermion play roles of up
and down quarks, and the composite particles (mesons) appear as triplets or
singlets of the $\mathrm{SU}(2)_{V}/\mathbb{Z}_{2}$ isospin symmetry. We used
the tensor network, in particular, DMRG for numerical simulation.

We obtained the masses of the pion, sigma, and eta meson by the three methods,
and the results are summarized in Fig. 10. We find that the results are
roughly consistent with each other taking into account possible systematic
errors for each method, such as the continuum and infinite-volume limits. In
addition, all the results show the relation $M_{\pi}<M_{\sigma}<M_{\eta}$,
which agrees with the analytic prediction by the bosonization technique. The
mass of the eta meson $M_{\eta}\sim 0.9$ is consistent with
$M_{\eta}=\mu+O(m)$ since $\mu\sim 0.8$ and $m=0.1$ in the current setup. We
also find that the masses of the pion and sigma mesons satisfy the WKB-based
formula (19), $M_{\sigma}/M_{\pi}=\sqrt{3}$, within not more than a 5%
deviation. It is, honestly, very surprising that the semiclassical analysis
gives the almost correct answer outside the range of its validity, and it
would be theoretically interesting to uncover the reason behind its success.

Figure 10:  The masses of the pion, sigma, and eta meson obtained by the three
independent methods are compared. Each result is obtained with the given
finite lattice spacing. We also put the error bar of the fitting error for the
correlation-function scheme, but it is too small to be seen.

Let us discuss the advantages and difficulties of each method and the
potential applications to other models. The first one, the correlation-
function scheme, is the straightforward generalization of the technique in
Lagrangian formalism. The advantage of this method is a wide range of
applicability to various models. We can obtain the meson masses from
correlation functions on a lattice with any dimensions, volume, and boundary
condition. Furthermore, the correlation function accepts the off-diagonal
element such as $\Braket{\mathcal{O}(x)\mathcal{O}^{\prime}(y)}$. This feature
will be useful when we turn on $\theta\neq 0$ in the 2-flavor Schwinger model.
The reason is that the meson operators become nontrivial mixtures of
$S_{f}(n)$ and $PS_{f}(n)$ depending on $\theta$. In this case, we need to
measure the correlation matrix of the operators and diagonalize it to extract
the mode of each meson. However, our numerical results suggest that the bond
dimension of MPS has to be sufficiently large to reproduce the correct
asymptotic behavior of the correlation function. In particular, the
computational cost increases rapidly as the system approaches a gapless phase,
for example, $m\sim 0$ or $\theta\sim\pi$. Thus, the tensor network (MPS) is
not an efficient approach to computing the mass spectrum by using correlation
functions.777 It is possible that other types of tensor networks, such as
MERA, may be still useful in this method. On the other hand, an ideal quantum
computer is free from such a restriction of the bound dimension. Thus, the
correlation function may be the first option in the era of practical quantum
computation of field theories in this sense, though to avoid the finite volume
effect for the two-point function we need a sizable scale computer.

The second method, the one-point-function scheme, makes good use of the
boundary effect rather than eliminating it. The results turn out to be
insensitive to the bond dimension, and thus we have to increase neither the
lattice size nor the bond dimension so much. Furthermore, the evaluation of
the local one-point function is generally easier than that of the long-range
correlation function. Thus, this is the most economical one among the three
methods. We note, however, that we have to prepare suitable boundary
conditions, such as defects, impurities, or external fields, to compute the
mass spectrum with this one-point-function scheme, which requires good
physical insights for the system of interest. In our case, the open boundary
at $\theta=0$ can be regarded as a source of the iso-singlet mesons, $\sigma$,
and $\eta$, but we have to set $\theta=2\pi$ to induce the boundary excitation
as a source of the iso-triplet mesons, $\pi_{a}$. We should also note that we
cannot obtain information on the off-diagonal correlators in the one-point-
function scheme. When $\theta=0$, the off-diagonal correlators are unimportant
because $\pi$, $\sigma$, and $\eta$ have different quantum numbers, but they
should become important at generic values of $\theta$ because the $G$ parity
is no longer a good quantum number.

The third method, the dispersion-relation scheme, is the distinctive strategy
of Hamiltonian formalism. We can obtain various states heuristically without
knowing what kind of mesons appear in the spectrum. We found the scattering
state unexpectedly by this method. Once we generate the excited states, it is
straightforward to measure various observables such as energy, momentum, and
quantum numbers. The states are identified by using these pieces of
information. Furthermore, we can investigate the wave function to distinguish
the $s$-wave or $p$-wave states. In this method, however, it is difficult to
increase the system size or the spatial dimensions. The reason is that we have
to generate an increasing number of states to search for different mesons. For
example, in our setup, we encounter the $3\times 4$ states of the pion before
obtaining the sigma meson at $k=13$. The momentum $K$ is discretized as $K\sim
2\pi\kappa/L$ for $\kappa=1,2,\cdots$ in the finite system with the size $L$.
If $L$ is increased, the number of pion states in a certain range of energy
grows up. Thus, we have to generate more excited states to reach the state of
the sigma meson. As for higher dimensions, there are momentum excitations in
each spatial direction, which result in an additional degeneracy. We expect
that there is a way to avoid this issue by modifying the strategy. For
example, if we are interested in a specific meson, it is more effective
generating excited states with a constraint on the quantum number to skip
mesons out of interest.

In this work, we have computed the mass spectrum at $\theta=0$. We note that
we have neglected many systematic errors, and thus there is plenty of room for
improvement. As a physics, extending our investigation to $\theta\neq 0$
should be interesting, where the sign problem arises in naive applications of
Monte Carlo simulations. The presence of $\theta$ introduces some differences
compared to the $\theta=0$ case. Firstly, the mass of the pion, which
corresponds to the gap of the system, decreases as $\theta\to\pi$.
Consequently, we may need to increase the bound dimension of MPS, leading to
higher computational costs. Secondly, the parity and $G$-parity are no longer
good quantum numbers for $\theta\not=0$, and the scalar and pseudo-scalar
operators have a nontrivial mixture. To handle this situation, we should
measure the correlation matrix between these operators and diagonalize it.
Although distinguishing the excited states, especially $\sigma$ and $\eta$,
seems to become tricky, exploring the changes in the spectrum promises
intriguing insights. Needless to say, it is very desirable that future
developments of these techniques eventually enable us to compute the hadron
spectrum of $4$d strongly-coupled gauge theories having the sign problem in
the conventional Monte Carlo methods.

###### Acknowledgements.

We would like to thank S. Aoki, M. Honda, T. Nishino, and K. Okunishi for
their useful discussions. The numerical calculations were carried out on XC40
at YITP in Kyoto University and the PC clusters at RIKEN iTHEMS. The work of
A. M. is supported by FY2022 Incentive Research Projects of RIKEN. The work of
E. I. is supported by JST PRESTO Grant Number JPMJPR2113, JST Grant Number
JPMJPF2221, JSPS KAKENHI (S) Grant number 23H05439, JSPS Grant-in-Aid for
Transformative Research Areas (A) JP21H05190, and Program for Promoting
Researches on the Supercomputer Fugaku” (Simulation for basic science:
approaching the new quantum era) Grant number JPMXP1020230411. The work of Y.
T. is supported by JSPS KAKENHI Grant number, 22H01218. This work is supported
by Center for Gravitational Physics and Quantum Information (CGPQI) at YITP.

##  Appendix A Operators in the spin representation

In this appendix, we show the spin representations of the Hamiltonian and
operators defined in Section 3 after the Jordan-Winger transformation (39) and
(40). For later convenience, we first show the transformation of some local
operators,

$\chi_{f,n}^{\dagger}\chi_{f,n}=\sigma_{f,n}^{+}\sigma_{f,n}^{-}=\frac{\sigma_{f,n}^{z}+1}{2},$
(89)
$\chi_{1,n}^{\dagger}\chi_{1,n+1}-\chi_{1,n+1}^{\dagger}\chi_{1,n}=\sigma_{1,n}^{+}\sigma_{2,n}^{z}\sigma_{1,n+1}^{-}-\sigma_{1,n}^{-}\sigma_{2,n}^{z}\sigma_{1,n+1}^{+},$
(90)
$\chi_{2,n}^{\dagger}\chi_{2,n+1}-\chi_{2,n+1}^{\dagger}\chi_{2,n}=\sigma_{2,n}^{+}\sigma_{1,n+1}^{z}\sigma_{2,n+1}^{-}-\sigma_{2,n}^{-}\sigma_{1,n+1}^{z}\sigma_{2,n+1}^{+}.$
(91)

The product of $\sigma^{z}$ in the Jordan-Winger transformation mostly cancels
each other in the fermion bilinears. Let us start with the Hamiltonian. Using
the relation above, the gauge part $H_{J}$ (33) is transformed as

$H_{J}=\frac{J}{4}\sum_{n=0}^{N-2}\left[\sum_{f=1}^{N_{f}}\sum_{k=0}^{n}\sigma_{f,k}^{z}+N_{f}\frac{(-1)^{n}+1}{2}+\frac{\theta}{\pi}\right]^{2}.$
(92)

The the fermion kinetic term $H_{w}$ (34) and the mass term $H_{m}$ (35) are
given by

$\displaystyle H_{w}$
$\displaystyle=-iw\sum_{n=0}^{N-2}\left(\sigma_{1,n}^{+}\sigma_{2,n}^{z}\sigma_{1,n+1}^{-}-\sigma_{1,n}^{-}\sigma_{2,n}^{z}\sigma_{1,n+1}^{+}\right.$
$\displaystyle\hphantom{=-iw\sum_{n=0}^{N-2}}\left.+\sigma_{2,n}^{+}\sigma_{1,n+1}^{z}\sigma_{2,n+1}^{-}-\sigma_{2,n}^{-}\sigma_{1,n+1}^{z}\sigma_{2,n+1}^{+}\right),$
(93)
$H_{m}=\frac{m_{\mathrm{lat}}}{2}\sum_{f=1}^{N_{f}}\sum_{n=0}^{N-1}(-1)^{n}\sigma_{f,n}^{z}+\frac{m_{\mathrm{lat}}}{2}N_{f}\frac{1-(-1)^{N}}{2}.$
(94)

Then the total Hamiltonian is a sum of them

$H=H_{J}+H_{w}+H_{m}.$ (95)

For practical use, $H_{J}$ can be decomposed into the quadratic term of
$\sigma^{z}$, the linear term of $\sigma^{z}$, and the constant term by
expanding the square. They can be summarized as follows:

$H_{J}=H_{J}^{(2)}+H_{J}^{(1)}+H_{J}^{(0)},$ (96)
$H_{J}^{(2)}=\frac{J}{2}\sum_{f=1}^{N_{f}}\sum_{j=0}^{N-3}\sum_{k=j+1}^{N-2}(N-k-1)\sigma_{f,j}^{z}\sigma_{f,k}^{z}+\frac{J}{4}\sum_{f\neq
f^{\prime}}\sum_{n=0}^{N-2}\sum_{j,k=0}^{n}\sigma_{f,j}^{z}\sigma_{f^{\prime},k}^{z},$
(97)
$H_{J}^{(1)}=\frac{J}{2}\sum_{f=1}^{N_{f}}\sum_{k=0}^{N-2}\left[\left(\frac{N_{f}}{2}+\frac{\theta}{\pi}\right)(N-k-1)+\frac{N_{f}}{2}\frac{(-1)^{N}+(-1)^{k}}{2}\right]\sigma_{f,k}^{z},$
(98) $\displaystyle H_{J}^{(0)}$
$\displaystyle=\frac{JN_{f}}{4}\frac{N(N-1)}{2}$
$\displaystyle+\frac{JN_{f}}{2}\left(\frac{N_{f}}{4}+\frac{\theta}{2\pi}\right)\left[\frac{(-1)^{N}-1}{2}+N\right]+J\left(\frac{\theta}{2\pi}\right)^{2}(N-1).$
(99)

The spin Hamiltonian contains the non-local interactions which come from the
Gauss law. It is not obvious whether the ground state can be described
efficiently by MPS.

Next, we map the observables by the Jordan-Winger transformation. The local
scalar condensate (48) and the pseudo-scalar condensate (49) are transformed
as

$S_{f}(n)=\frac{1}{8a}(-1)^{n}(-\sigma_{f,n-1}^{z}+2\sigma_{f,n}^{z}-\sigma_{f,n+1}^{z}),$
(100) $\displaystyle PS_{1}(n)=\frac{i}{4a}(-1)^{n}$
$\displaystyle\left(\sigma_{1,n-1}^{+}\sigma_{2,n-1}^{z}\sigma_{1,n}^{-}-\sigma_{1,n-1}^{-}\sigma_{2,n-1}^{z}\sigma_{1,n}^{+}\right.$
$\displaystyle\left.-\sigma_{1,n}^{+}\sigma_{2,n}^{z}\sigma_{1,n+1}^{-}+\sigma_{1,n}^{-}\sigma_{2,n}^{z}\sigma_{1,n+1}^{+}\right),$
(101) $\displaystyle PS_{2}(n)=\frac{i}{4a}(-1)^{n}$
$\displaystyle\left(\sigma_{2,n-1}^{+}\sigma_{1,n}^{z}\sigma_{2,n}^{-}-\sigma_{2,n-1}^{-}\sigma_{1,n}^{z}\sigma_{2,n}^{+}\right.$
$\displaystyle\left.-\sigma_{2,n}^{+}\sigma_{1,n+1}^{z}\sigma_{2,n+1}^{-}+\sigma_{2,n}^{-}\sigma_{1,n+1}^{z}\sigma_{2,n+1}^{+}\right).$
(102)

We can also map the isospin operators (50), (52), and (53) as follows:

$J_{z}=\frac{1}{4}\sum_{n=0}^{N-1}(\sigma_{1,n}^{z}-\sigma_{2,n}^{z}),$ (103)
$J_{+}=i\sum_{n=0}^{N-1}\sigma_{1,n}^{+}\sigma_{2,n}^{-},$ (104)
$J_{-}=-i\sum_{n=0}^{N-1}\sigma_{2,n}^{+}\sigma_{1,n}^{-}.$ (105)

Then it is straightforward to construct these operators as MPO.

Finally, we consider the Jordan-Winger transformation of the total momentum
operator (65). Each term in the sum is transformed as follows:

$\displaystyle\chi_{1,n-1}^{\dagger}\chi_{1,n+1}-\chi_{1,n+1}^{\dagger}\chi_{1,n-1}$
$\displaystyle=-\sigma_{1,n-1}^{+}\sigma_{2,n-1}^{z}\sigma_{1,n}^{z}\sigma_{2,n}^{z}\sigma_{1,n+1}^{-}+\sigma_{1,n-1}^{-}\sigma_{2,n-1}^{z}\sigma_{1,n}^{z}\sigma_{2,n}^{z}\sigma_{1,n+1}^{+},$
(106)
$\displaystyle\chi_{2,n-1}^{\dagger}\chi_{2,n+1}-\chi_{2,n+1}^{\dagger}\chi_{2,n-1}$
$\displaystyle=-\sigma_{2,n-1}^{+}\sigma_{1,n}^{z}\sigma_{2,n}^{z}\sigma_{1,n+1}^{z}\sigma_{2,n+1}^{-}+\sigma_{2,n-1}^{-}\sigma_{1,n}^{z}\sigma_{2,n}^{z}\sigma_{1,n+1}^{z}\sigma_{2,n+1}^{+}.$
(107)

Thus, the total momentum is given by the combination of five Pauli matrices,

$\displaystyle K=\frac{i}{4a}\sum_{n=1}^{N-2}$
$\displaystyle\left(\sigma_{1,n-1}^{-}\sigma_{2,n-1}^{z}\sigma_{1,n}^{z}\sigma_{2,n}^{z}\sigma_{1,n+1}^{+}-\sigma_{1,n-1}^{+}\sigma_{2,n-1}^{z}\sigma_{1,n}^{z}\sigma_{2,n}^{z}\sigma_{1,n+1}^{-}\right.$
$\displaystyle\left.+\sigma_{2,n-1}^{-}\sigma_{1,n}^{z}\sigma_{2,n}^{z}\sigma_{1,n+1}^{z}\sigma_{2,n+1}^{+}-\sigma_{2,n-1}^{+}\sigma_{1,n}^{z}\sigma_{2,n}^{z}\sigma_{1,n+1}^{z}\sigma_{2,n+1}^{-}\right).$
(108)

##  Appendix B Charge conjugation operator in the 1-flavor Schwinger model

The charge conjugation operator $C$ defined by (58) does not commute with the
Hamiltonian with the open boundary condition. This is because the charge
conjugation for the staggered fermion must incorporate the one-unit lattice
translation, and thus it is not an on-site symmetry in our regularization
scheme. As a consequence, the expectation value of $C$ does not become $\pm
1$. However, this is an important quantum number to diagnose the type of
mesons, and we have assumed in Section 5.3 that we can diagnose the quantum
number by the sign of $\langle C\rangle$. Although we have no theoretical
justifications for this prescription, let us test it in the $1$-flavor
Schwinger model to give some evidence for its reasonableness.

Figure 11:  The expectation values of the charge conjugation $C$ for the
ground state (left) and the first excited state (right) are plotted against
the lattice spacing $a$. The number of lattice sites $N$ is chosen to fix the
physical length $L=(N-1)a$. Each symbol corresponds to a different value of
$L$. We set $\theta=0$ and $m=0.125$ in this analysis. The fitting results by
the quadratic are also shown by the solid lines.

First, we investigate the behavior of $C$ in the continuum limit. We generated
the MPS of the ground state and the 1st excited state of the 1-flavor
Schwinger model at $\theta=0$ by DMRG. The lattice spacing $a$ is changed
around $0.1\lesssim a\lesssim 0.25$. The number of lattice sites $N$ is chosen
to fix the physical system size $L=(N-1)a$. We compute the expectation values
of $C$ for these MPS. The results are shown in Fig. 11. The different symbols
correspond to the results for different $L$ in the plots. We fitted the data
points for each $L$ by the quadratic function $f(a)=c_{0}+c_{1}a+c_{2}a^{2}$.
The fitting results are also plotted in Fig. 11 by the solid lines. For
$L=49.8$, we obtained the continuum limit $\Braket{C}_{a\rightarrow
0}=0.321(3)$ for the ground state and $-0.320(3)$ for the 1st excited state.
The results with the other $L$ agree with these values within the error. Thus,
we confirmed that the expectation value of $C$ is a finite value in the
continuum limit and is not sensitive to $L$.

Next, let us discuss the effect of the boundary on $C$. We consider a further
simplified model, the free fermion on the periodic lattice. The model is
obtained from the 1-flavor Schwinger model with the periodic boundary
condition by setting $g=0$ and adding the hopping term between $n=0$ and
$n=N-1$ site. In fact, it is hard to adopt the p.b.c. in the current DMRG
method due to the artificial long-range interaction between both ends of MPS.
Thus, we choose small sizes of the lattice $N=20$ and $40$ for this analysis.
The corresponding lattice spacings are $a=0.2$ and $0.1$ for the fixed
physical length $L=Na=4$. We generate the ground state and the excited states
up to the level $k=4$. The four excited states turned out to be degenerated.
Thus, we compute $\Braket{C}_{k,k^{\prime}}$ including the off-diagonal
elements, and diagonalize the result as the $4\times 4$ matrix. The
eigenvalues are shown in Fig. 12. We found that $\Braket{C}=1$ for the ground
state and $\Braket{C}=\pm 1,\alpha\pm i\beta$ for the excited states. These
complex values satisfy $|\alpha|^{2}+|\beta|^{2}=1$ as we can see in the plot.
The imaginary part $\beta$ becomes smaller as $a$ is decreased, which suggests
that we will obtain $\Braket{C}\rightarrow\pm 1$ in the continuum limit.

Figure 12:  The expectation values of $C$ are plotted on the complex plane.
The different symbols represent the results for the different values of the
spacing $a$. All the data points turned out to be on the unit circle.

##  Appendix C Arrangement of flavors on MPS

In the spin representation of the Hamiltonian (95) of the 2-flavor Schwinger
model, each spin has the site index $n$ and the flavor index $f$. To apply
DMRG, we arrange these spins on the 1d lattice with the single index
$(f,n)\rightarrow i$. Although the ordering of the indices does not affect the
physics, it can affect the necessary bound dimensions, and thus calculation
cost depends on it. In this work, we assign the index $i$ to $(f,n)$ as

$i=N_{f}n+f-1=0,1,\cdots,N_{f}N,$ (109)

which we call the staggered order in this section. In this arrangement,
different flavors at the same physical site are put closely with each other,
and this is important to control the bond dimension in the computation of
DMRG. Let us consider another choice for comparison,

$i=n+N(f-1),$ (110)

which we name the flavor order here. In this case, we first arrange one of the
flavors and then start to arrange the next one, so the flavor degrees of
freedom at the same physical sites are separated by $N$, and this clearly
violates the above important criterion.

Figure 13:  The effective bond dimension $D_{\mathrm{eff}}$ is plotted against
the number of sweeps $N_{\mathrm{sweep}}$ for the flavor order (left) and the
staggered order (right). The vertical axis of the left panel is in log scale,
whereas the axis of the right panel is in linear scale. The lattice spacing
and the fermion mass are set to $a=0.2$ and $m=0.1$. Figure 14:  The effective
bond dimension $D_{\mathrm{eff}}$ after 20 sweeps is plotted against the
system size $N$ in log-log scale. The result grows exponentially with $N$ for
the flavor order whereas it is saturated for the staggered order. The lattice
spacing and the fermion mass are set to $a=0.2$ and $m=0.1$.

In these two cases, we compare the efficiency of the MPS to represent the
ground state in the gapped phase $\theta=0$. We obtain the ground state by
DMRG and investigate the largest bond dimension in the MPS, called the
effective bond dimension $D_{\mathrm{eff}}$. The results are plotted against
the number of sweeps $N_{\mathrm{sweep}}$ in Fig. 13 for various lattice sizes
$N$. We found that $D_{\mathrm{eff}}$ converges around $O(10)$ sweeps for both
cases. However, the dependence on $N$ is totally different. For the flavor
order, the final value of $D_{\mathrm{eff}}$ increases exponentially with $N$,
which is caused by artificial long-range interaction between the two flavors.
On the other hand for the staggered order, the final value is saturated for
sufficiently large $N$. To show these behaviors, we plot the final values of
$D_{\mathrm{eff}}$ after 20 sweeps against $N$ in Fig. 14.

According to Fig. 14, the bond dimension seems to saturate in the case of the
staggered order as $N\to\infty$. Since the $\ln D_{\mathrm{eff}}$ gives the
upper bound for the entanglement entropy, this constant behavior is expected
to be the optimal one for the $1+1$d gapped systems. On the other hand,
$D_{\mathrm{eff}}$ grows exponentially fast for the flavor ordered as
$N\to\infty$. We suspect that this is because the flavor order puts the
entangled flavors in separate locations. If we cut the system into two pieces
in terms of $i$ with the flavor ordering, the $O(N)$ entangled pairs are cut,
and thus the entanglement entropy becomes $O(N)$, which is consistent with the
exponential behavior of $D_{\mathrm{eff}}$.

Therefore, we adopt the staggered order in the whole analysis of this work.

##  Appendix D Correlation function in the 1-flavor Schwinger model

We test the validity of the correlation-function scheme in Section 5.1 by
examining the correlation function in the 1-flavor Schwinger model. When the
fermion is massless $m=0$, the model can be analytically solvable and it is
equivalent to the free massive boson with mass $\mu^{\prime}=g/\sqrt{\pi}$.
Thus, this is a good benchmark and we compare the numerical result of DMRG
with the analytical answer.

As an analogy of the pseudo scalar meson in the 2-flavor Schwinger model, we
consider the pseudo-scalar operator $PS=-i\bar{\psi}\gamma^{5}\psi$. The
results of the correlation function $\Braket{PS(x)PS(y)}$ are shown in the
left panel of Fig. 15. Here, the data with different colors are obtained with
the different values of the cutoff parameter $\varepsilon$. The corresponding
effective masses (3-point average) are also plotted in the right panel of Fig.
15, where we can see the significant $\varepsilon$ dependence.

Figure 15:  (Left) The correlation function $\ln\Braket{PS(x)PS(y)}$ is
plotted against the distance $r=|x-y|$ for various values of $\varepsilon$.
The number of lattice sites is $N=400$ and the lattice spacing $a$ is
determined so that $L=a(N-1)=79.8$. (Right) The effective mass
$M_{\mathrm{eff}}(r)$ (3-point average) calculated from the correlation
function in the left panel is plotted against $r$.

To see the $1/r$ correction of the effective mass, we plot
$M_{\mathrm{eff}}(r)$ against $1/r$ in Fig. 16. Then we found that the result
approaches the expected asymptotic behavior
$M_{\mathrm{eff}}(r)\sim\alpha/r+M$ only if the cutoff $\varepsilon$ is
sufficiently small. We fitted the data points for $\varepsilon=10^{-16}$ by
$\alpha/r+M$ in the range $0.06\leq 1/r\leq 0.2$ and obtained $M=0.5677(5)$
and $\alpha=0.446(4)$. Here the systematic error from the uncertainty of the
fitting range is evaluated as explained in Section 5.1. We note that this is
the result on the finite lattice before taking the continuum limit, but it
turned out to be close to the exact value $M=g/\sqrt{\pi}\approx 0.56419$ of
the continuum theory.

Therefore, it is quite important to discuss the cutoff (or bond-dimension)
dependence especially when we use the correlation-function scheme. Indeed, if
we naively read the plateau value at $\varepsilon=10^{-10}$, we got an
incorrect answer $M\sim 0.63$ without observing the $1/\sqrt{r}$ contribution
in the Yukawa potential at all.

Figure 16:  The effective mass $M_{\mathrm{eff}}(r)$ is plotted against $1/r$.
The data points for $\varepsilon=10^{-16}$ are fitted by $\alpha/r+M$ inside
the region $0.06\leq 1/r\leq 0.2$. The fitting result is depicted by the
shaded band with systematic error. The exact mass of the pseudo scalar
$g/\sqrt{\pi}$ is also shown by the horizontal broken line.

## References

  * (1) Flavour Lattice Averaging Group (FLAG) collaboration, _FLAG Review 2021_ , _Eur. Phys. J. C_ 82 (2022) 869 [2111.09849].
  * (2) S. Borsanyi, Z. Fodor, C. Hoelbling, S.D. Katz, S. Krieg and K.K. Szabo, _Full result for the QCD equation of state with 2+1 flavors_ , _Phys. Lett. B_ 730 (2014) 99 [1309.5258].
  * (3) HotQCD collaboration, _Equation of state in ( 2+1 )-flavor QCD_ , _Phys. Rev. D_ 90 (2014) 094503 [1407.6387].
  * (4) N. Shibata, K. Ueda, T. Nishino and C. Ishii, _Friedel oscillations in the one-dimensional kondo lattice model_ , _Phys. Rev. B_ 54 (1996) 13495.
  * (5) N. Shibata, K. Ueda, T. Nishino and C. Ishii, _Large fermi surface of the one-dimensional kondo lattice model observed by friedel oscillations_ , _Physica B: Condensed Matter_ 230-232 (1997) 1024.
  * (6) B. Pirvu, J. Haegeman and F. Verstraete, _Matrix product state based algorithm for determining dispersion relations of quantum spin chains with periodic boundary conditions_ , _Physical Review B_ 85 (2012) .
  * (7) J. Haegeman, B. Pirvu, D.J. Weir, J.I. Cirac, T.J. Osborne, H. Verschelde et al., _Variational matrix product ansatz for dispersion relations_ , _Physical Review B_ 85 (2012) .
  * (8) J. Haegeman, S. Michalakis, B. Nachtergaele, T.J. Osborne, N. Schuch and F. Verstraete, _Elementary excitations in gapped quantum spin systems_ , _Physical Review Letters_ 111 (2013) .
  * (9) J.S. Schwinger, _Gauge Invariance and Mass. 2._ , _Phys. Rev._ 128 (1962) 2425.
  * (10) J.H. Lowenstein and J.A. Swieca, _Quantum electrodynamics in two-dimensions_ , _Annals Phys._ 68 (1971) 172.
  * (11) A. Casher, J.B. Kogut and L. Susskind, _Vacuum polarization and the absence of free quarks_ , _Phys. Rev. D_ 10 (1974) 732.
  * (12) S.R. Coleman, R. Jackiw and L. Susskind, _Charge Shielding and Quark Confinement in the Massive Schwinger Model_ , _Annals Phys._ 93 (1975) 267.
  * (13) S.R. Coleman, _More About the Massive Schwinger Model_ , _Annals Phys._ 101 (1976) 239.
  * (14) N.S. Manton, _The Schwinger Model and Its Axial Anomaly_ , _Annals Phys._ 159 (1985) 220.
  * (15) J.E. Hetrick and Y. Hosotani, _QED ON A CIRCLE_ , _Phys. Rev._ D38 (1988) 2621.
  * (16) C. Jayewardena, _SCHWINGER MODEL ON S(2)_ , _Helv. Phys. Acta_ 61 (1988) 636.
  * (17) I. Sachs and A. Wipf, _Finite temperature Schwinger model_ , _Helv. Phys. Acta_ 65 (1992) 652 [1005.1822].
  * (18) C. Adam, _Instantons and vacuum expectation values in the Schwinger model_ , _Z. Phys._ C63 (1994) 169.
  * (19) C. Adam, _The Dyson-Schwinger equations in the instanton vacuum of the Schwinger model_ , _Czech. J. Phys._ 46 (1996) 893 [hep-ph/9501273].
  * (20) J.E. Hetrick, Y. Hosotani and S. Iso, _The Massive multi - flavor Schwinger model_ , _Phys. Lett._ B350 (1995) 92 [hep-th/9502113].
  * (21) R. Narayanan, _QED at a finite chemical potential_ , _Phys. Rev._ D86 (2012) 087701 [1206.1489].
  * (22) R. Narayanan, _Two flavor massless Schwinger model on a torus at a finite chemical potential_ , _Phys. Rev._ D86 (2012) 125008 [1210.3072].
  * (23) R. Lohmayer and R. Narayanan, _Phase structure of two-dimensional QED at zero temperature with flavor-dependent chemical potentials and the role of multidimensional theta functions_ , _Phys. Rev._ D88 (2013) 105030 [1307.4969].
  * (24) Y. Tanizaki and M. Tachibana, _Multi-flavor massless QED 2 at finite densities via Lefschetz thimbles_, _JHEP_ 02 (2017) 081 [1612.06529].
  * (25) M.C. Bañuls, K. Cichy, K. Jansen and J.I. Cirac, _The mass spectrum of the Schwinger model with Matrix Product States_ , _JHEP_ 11 (2013) 158 [1305.3765].
  * (26) M.C. Bañuls, K. Cichy, J.I. Cirac, K. Jansen and H. Saito, _Thermal evolution of the Schwinger model with Matrix Product Operators_ , _Phys. Rev. D_ 92 (2015) 034519 [1505.00279].
  * (27) M.C. Bañuls, K. Cichy, K. Jansen and H. Saito, _Chiral condensate in the Schwinger model with Matrix Product Operators_ , _Phys. Rev. D_ 93 (2016) 094512 [1603.05002].
  * (28) B. Buyens, F. Verstraete and K. Van Acoleyen, _Hamiltonian simulation of the Schwinger model at finite temperature_ , _Phys. Rev. D_ 94 (2016) 085018 [1606.03385].
  * (29) 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 (2017) 114501 [1612.00739].
  * (30) M.C. Bañuls, K. Cichy, J.I. Cirac, K. Jansen and S. Kühn, _Density Induced Phase Transitions in the Schwinger Model: A Study with Matrix Product States_ , _Phys. Rev. Lett._ 118 (2017) 071601 [1611.00705].
  * (31) L. Funcke, K. Jansen and S. Kühn, _Topological vacuum structure of the Schwinger model with matrix product states_ , _Phys. Rev. D_ 101 (2020) 054507 [1908.00551].
  * (32) B. Chakraborty, M. Honda, T. Izubuchi, Y. Kikuchi and A. Tomiya, _Classically emulated digital quantum simulation of the Schwinger model with a topological term via adiabatic state preparation_ , _Phys. Rev. D_ 105 (2022) 094503 [2001.00485].
  * (33) 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 (2022) 014504 [2105.03276].
  * (34) M. Honda, E. Itou, Y. Kikuchi and Y. Tanizaki, _Negative string tension of a higher-charge Schwinger model via digital quantum simulation_ , _PTEP_ 2022 (2022) 033B01 [2110.14105].
  * (35) M. Honda, E. Itou and Y. Tanizaki, _DMRG study of the higher-charge Schwinger model and its ’t Hooft anomaly_ , _JHEP_ 11 (2022) 141 [2210.04237].
  * (36) A. Tomiya, _Schwinger model at finite temperature and density with beta VQE_ , 2205.08860.
  * (37) L. Funcke, K. Jansen and S. Kühn, _Exploring the CP-Violating Dashen Phase in the Schwinger Model with Tensor Networks_ , 2303.03799.
  * (38) R. Dempsey, I.R. Klebanov, S.S. Pufu, B.T. Søgaard and B. Zan, _Phase Diagram of the Two-Flavor Schwinger Model at Zero Temperature_ , 2305.04437.
  * (39) D.E. Kharzeev and Y. Kikuchi, _Real-time chiral dynamics from a digital quantum simulation_ , _Phys. Rev. Res._ 2 (2020) 023342 [2001.00698].
  * (40) W.A. de Jong, K. Lee, J. Mulligan, M. Płoskoń, F. Ringer and X. Yao, _Quantum simulation of nonequilibrium dynamics and thermalization in the Schwinger model_ , _Phys. Rev. D_ 106 (2022) 054508 [2106.08394].
  * (41) N.H. Nguyen, M.C. Tran, Y. Zhu, A.M. Green, C.H. Alderete, Z. Davoudi et al., _Digital Quantum Simulation of the Schwinger Model and Symmetry Protection with Trapped Ions_ , _PRX Quantum_ 3 (2022) 020324 [2112.14262].
  * (42) L. Nagano, A. Bapat and C.W. Bauer, _Quench dynamics of the Schwinger model via variational quantum algorithms_ , 2302.10933.
  * (43) H. Fukaya and T. Onogi, _Lattice study of the massive Schwinger model with theta term under Luscher’s ’admissibility’ condition_ , _Phys. Rev. D_ 68 (2003) 074503 [hep-lat/0305004].
  * (44) C. Gattringer, T. Kloiber and V. Sazonov, _Solving the sign problems of the massless lattice Schwinger model with a dual formulation_ , _Nucl. Phys. B_ 897 (2015) 732 [1502.05479].
  * (45) C. Gattringer, T. Kloiber and M. Müller-Preussker, _Dual simulation of the two-dimensional lattice U(1) gauge-Higgs model with a topological term_ , _Phys. Rev._ D92 (2015) 114508 [1508.00681].
  * (46) C. Gattringer, D. Göschl and T. Sulejmanpasic, _Dual simulation of the 2d U(1) gauge Higgs model at topological angle $\theta=\pi\,$: Critical endpoint behavior_, _Nucl. Phys._ B935 (2018) 344 [1807.07793].
  * (47) K. Harada, T. Sugihara, M.-a. Taniguchi and M. Yahiro, _The Massive Schwinger model with SU(2)-f on the light cone_ , _Phys. Rev. D_ 49 (1994) 4226 [hep-th/9309128].
  * (48) M. Fishman, S.R. White and E.M. Stoudenmire, _The ITensor Software Library for Tensor Network Calculations_ , _SciPost Phys. Codebases_ (2022) 4.
  * (49) T.D. Lee and C.N. Yang, _Charge conjugation, a new quantum numberg, and selection rules concerning a nucleon-antinucleon system_ , _Il Nuovo Cimento (1955-1965)_ 3 (1956) 749.
  * (50) J. Frohlich and E. Seiler, _The Massive Thirring-Schwinger Model (QED in Two-Dimensions): Convergence of Perturbation Theory and Particle Structure_ , _Helv. Phys. Acta_ 49 (1976) 889.
  * (51) X. Chen, Z.-C. Gu and X.-G. Wen, _Classification of gapped symmetric phases in one-dimensional spin systems_ , _Phys. Rev. B_ 83 (2011) 035107 [1008.3745].
  * (52) A. Kapustin, _Bosonic Topological Insulators and Paramagnets: a view from cobordisms_ , 1404.6659.
  * (53) A. Kapustin, _Symmetry Protected Topological Phases, Anomalies, and Cobordisms: Beyond Group Cohomology_ , 1403.1467.
  * (54) T. Misumi, Y. Tanizaki and M. Ünsal, _Fractional $\theta$ angle, ’t Hooft anomaly, and quantum instantons in charge-$q$ multi-flavor Schwinger model_, _JHEP_ 07 (2019) 018 [1905.05781].
  * (55) F.D.M. Haldane, _Nonlinear field theory of large spin Heisenberg antiferromagnets. Semiclassically quantized solitons of the one-dimensional easy Axis Neel state_ , _Phys. Rev. Lett._ 50 (1983) 1153.
  * (56) I. Affleck and E.H. Lieb, _A Proof of Part of Haldane’s Conjecture on Spin Chains_ , _Lett. Math. Phys._ 12 (1986) 57.
  * (57) F.D.M. Haldane, _O (3) Nonlinear sigma Model and the Topological Distinction between Integer- and Half-Integer-Spin Antiferromagnets in Two Dimensions_ , _Phys. Rev. Lett._ 61 (1988) 1029.
  * (58) I. Affleck, T. Kennedy, E.H. Lieb and H. Tasaki, _Rigorous Results on Valence Bond Ground States in Antiferromagnets_ , _Phys. Rev. Lett._ 59 (1987) 799.
  * (59) Z. Komargodski, A. Sharon, R. Thorngren and X. Zhou, _Comments on Abelian Higgs Models and Persistent Order_ , _SciPost Phys._ 6 (2019) 003 [1705.04786].
  * (60) Z. Komargodski, T. Sulejmanpasic and M. Unsal, _Walls, anomalies, and deconfinement in quantum antiferromagnets_ , _Phys. Rev._ B97 (2018) 054418 [1706.05731].
  * (61) M. Lajkó, K. Wamer, F. Mila and I. Affleck, _Generalization of the Haldane conjecture to SU(3) chains_ , _Nucl. Phys._ B924 (2017) 508 [1706.06598].
  * (62) Y. Tanizaki and T. Sulejmanpasic, _Anomaly and global inconsistency matching: $\theta$-angles, $SU(3)/U(1)^{2}$ nonlinear sigma model, $SU(3)$ chains and its generalizations_, _Phys. Rev._ B98 (2018) 115126 [1805.11423].
  * (63) S.R. Coleman, _The Quantum Sine-Gordon Equation as the Massive Thirring Model_ , _Phys. Rev. D_ 11 (1975) 2088.
  * (64) R.F. Dashen, B. Hasslacher and A. Neveu, _The Particle Spectrum in Model Field Theories from Semiclassical Functional Integral Techniques_ , _Phys. Rev. D_ 11 (1975) 3424.
  * (65) J.B. Kogut and L. Susskind, _Hamiltonian Formulation of Wilson’s Lattice Gauge Theories_ , _Phys. Rev. D_ 11 (1975) 395.
  * (66) L. Susskind, _Lattice Fermions_ , _Phys. Rev. D_ 16 (1977) 3031.
  * (67) 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 (2022) 043133 [2206.05308].
  * (68) S.R. White, _Density matrix formulation for quantum renormalization groups_ , _Phys. Rev. Lett._ 69 (1992) 2863.
  * (69) S.R. White, _Density-matrix algorithms for quantum renormalization groups_ , _Phys. Rev. B_ 48 (1993) 10345.
  * (70) U. Schollwöck, _The density-matrix renormalization group_ , _Reviews of Modern Physics_ 77 (2005) 259.
  * (71) U. Schollwöck, _The density-matrix renormalization group in the age of matrix product states_ , _Annals of Physics_ 326 (2011) 96.
  * (72) E. Stoudenmire and S.R. White, _Studying two-dimensional systems with the density matrix renormalization group_ , _Annual Review of Condensed Matter Physics_ 3 (2012) 111.
  * (73) M.L. Wall and L.D. Carr, _Out-of-equilibrium dynamics with matrix product states_ , _New Journal of Physics_ 14 (2012) 125015.