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

# DMRG study of the theta-dependent mass spectrum in the $2$-flavor Schwinger
model

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 study the $\theta$-dependent mass spectrum of the massive $2$-flavor
Schwinger model in the Hamiltonian formalism using the density-matrix
renormalization group(DMRG). The masses of the composite particles, the pion
and sigma meson, are computed by two independent methods. One is the improved
one-point-function scheme, where we measure the local meson operator coupled
to the boundary state and extract the mass from its exponential decay. Since
the $\theta$ term causes a nontrivial operator mixing, we unravel it by
diagonalizing the correlation matrix to define the meson operator. The other
is the dispersion-relation scheme, a heuristic approach specific to
Hamiltonian formalism. We obtain the dispersion relation directly by measuring
the energy and momentum of the excited states. The sign problem is
circumvented in these methods, and their results agree with each other even
for large $\theta$. We reveal that the $\theta$-dependence of the pion mass at
$m/g=0.1$ is consistent with the prediction by the bosonized model. We also
find that the mass of the sigma meson satisfies the semi-classical formula,
$M_{\sigma}/M_{\pi}=\sqrt{3}$, for almost all region of $\theta$. While the
sigma meson is a stable particle thanks to this relation, the eta meson is no
longer protected by the $G$-parity and becomes unstable for $\theta\neq 0$.

††preprint: YITP-24-79, RIKEN-iTHEMS-Report-24

##  1 Introduction

Lattice gauge theories give the foundation of strongly coupled quantum field
theories (QFTs). The Monte Carlo simulation based on its Euclidean path
integral has uncovered various non-perturbative aspects, including properties
of hadrons in quantum chromodynamics (QCD), and those results are successfully
comparable to experimental results FlavourLatticeAveragingGroupFLAG:2021npn .
One of the notable achievements in the Lattice Monte Carlo study is the
precise prediction of the hadron mass spectra with only a few input
parameters. Indeed, the mass spectra of the various hadrons cannot be
explained by the sum of the renormalized quark masses of the constituents, and
the contributions from interactions need to be evaluated non-perturbatively.
However, it is difficult to perform the Monte Carlo simulation when it suffers
from the sign problem, which typically appears in the non-zero chemical
potential or with the topological $\theta$ term in QCD Nagata:2021ugx . We
need to expand the scope of numerical studies of QFTs to uncover the new
physics hidden in such regions.

Quantum computing and tensor network methods based on the Hamiltonian
formalism of lattice gauge theories Kogut:1974ag  have recently gotten
attention as an alternative numerical approach. These methods directly
approximate wave functionals of low-energy states, so they are based on a
completely different criterion from that of the importance sampling. The
development of these complementary frameworks may enable us to obtain new
information that is difficult to extract from the Euclidean path integral.

With the current of the times, we have developed three distinct methods for
computing mass spectra of lattice Hamiltonian gauge theories in ref.
Itou:2023img : (1) The correlation-function scheme extracts the particle mass
by measuring the spatial correlator, which is similar to the conventional
Euclidean method. (2) The one-point-function scheme utilizes the boundary
effect as a source of the composite particles, and the mass is computed from
the exponential decay of the one-point function. (3) The dispersion-relation
scheme directly computes the dispersion relation by generating the low-lying
states at finite volume, where the states of each meson with the momentum
excitation can be identified via the quantum numbers. The last scheme, in
particular, is specific to the Hamiltonian formalism, unlike the other two. In
ref. Itou:2023img , we applied these methods to the $2$-flavor massive
Shcwinger model at $\theta=0$ using the density-matrix renormalization group
(DRMG) White:1992zz ; White:1993zza ; Schollw_ck_2005 ; Schollw_ck_2011 , and
we confirmed that the mass spectra obtained by these methods are almost
consistent with each other and also with theoretical predictions.

The Schwinger model Schwinger:1962tp  is the $(1+1)$d quantum electrodynamics,
showing charge confinement as 4d QCD. Its Euclidean formulation suffers from
the sign problem at nonzero $\theta$ angles, but many of its nonperturbative
properties are analytically calculable in the massless limit Lowenstein:1971fc
; Casher:1974vf ; Coleman:1975pw ; Manton:1985jm ; Hetrick:1988yg ;
Jayewardena:1988td ; Sachs:1991en ; Adam:1993fc ; Hetrick:1995wq ;
Narayanan:2012du ; Narayanan:2012qf ; Lohmayer:2013eka ; Tanizaki:2016xcu .
Numerically, the $\theta$-dependence of the pion mass for the $2$-flavor case
is computed using the reweighting technique for the conventional importance
sampling method, which gives reliable results up to $\theta/2\pi\lesssim 0.25$
and confirms some analytic predictions Fukaya:2003ph . However, the sign
problem becomes severe beyond that point and we need further different
approaches to numerically study the physics at larger $\theta$. Hence, the
Schwinger model provides an interesting testing ground for new computational
methods, and indeed there have been many studies on the model with nonzero
$\theta$ by the tensor network or quantum algorithm Banuls:2013jaa ;
Banuls:2015sta ; Banuls:2016lkq ; Buyens:2015tea ; Buyens:2016ecr ;
Buyens:2016hhu ; Buyens:2017crb ; Funcke:2019zna ; Chakraborty:2020uhf ;
Kharzeev:2020kgc ; Honda:2021aum ; Thompson:2021eze ; Honda:2021ovk ;
Honda:2022edn ; Nagano:2023uaq ; Funcke:2023lli ; Angelides:2023bme ;
Ikeda:2023zil ; Dempsey:2023gib ; Pedersen:2023asd ; Schmoll:2023eez ;
Angelides:2023noe ; Ghim:2024pxe ; Popov:2024ysn ; Kaikov:2024acw .111 As
another Hamiltonian-based approach, there is a numerical study with the light-
cone quantization Harada:1993va . The tensor network in the Lagrangian
formalism also provides a way to be liberated from importance sampling
Shimizu:2014uva ; Shimizu:2014fsa ; Butt:2019uul ; Yosprakob:2023tyr ; Az-
zahra:2024gqr . For a specific class of $(1+1)$d gauge theories including the
Schwinger model, the dual variable formulation Gattringer:2015nea ;
Gattringer:2015baa ; Gattringer:2018dlw ; Sulejmanpasic:2019ytl ;
Sulejmanpasic:2020lyq  or the bosonization Ohata:2023sqc ; Ohata:2023gru  can
be used to eliminate the sign problem in the Monte Carlo methods. Among these
options, we focus on the Hamiltonian formalism written by the spin variables
to fit the quantum computations and tensor networks.

In this paper, we extend our previous study Itou:2023img  to investigate the
physics of the $2$-flavor Schwinger model in the $\theta\neq 0$ regime. At
$\theta=0$, the mesons can be characterized by their quantum numbers $J^{PG}$,
where $J$, $P$, and $G$ are the isospin, parity, and $G$-parity, respectively.
There are three types of stable mesons; pion $\pi_{a}$ ($J^{PG}=1^{-+}$),
sigma meson $\sigma$ ($J^{PG}=0^{++}$), and eta meson $\eta$
($J^{PG}=0^{--}$). The low-energy spectrum has been studied analytically using
the bosonization technique Coleman:1976uz . Furthermore, the WKB-type
approximation predicts that the sigma-meson mass has the specific ratio with
the pion mass, $M_{\sigma}/M_{\pi}=\sqrt{3}$, and these are confirmed in our
previous work Itou:2023img . On the other hand, at $\theta\neq 0$, the
situation will be changed: First, the parity and $G$-parity are no longer good
quantum numbers at $\theta\neq 0$, and the symmetries can no longer
distinguish $\sigma$ and $\eta$. This causes the $\sigma$-$\eta$ mixing and
the decay channel $\eta\to\pi\pi$ also opens. Second, the $\theta$-dependence
of the pion mass is predicted to behave as $M_{\pi}(\theta)\sim
g|(m/g)\cos(\theta/2)|^{2/3}$ when the fermion mass is small Coleman:1976uz ,
and the mass of the sigma meson is also predicted to satisfy the relation
$M_{\sigma}(\theta)/M_{\pi}(\theta)=\sqrt{3}$, but it has not yet been
confirmed by the first-principle calculations because of the sign problem.
Moreover, substituting $\theta=\pi$ into these formulas, the spectrum is
expected to become massless. Indeed, it is suggested that the massive
$2$-flavor Schwinger model at $\theta=\pi$ shows almost conformal behavior and
is described by the $SU(2)_{1}$ Wess-Zumino-Witten (WZW) model with the
marginally relevant $J_{L}J_{R}$ deformation Coleman:1976uz ; Dempsey:2023gib
. From these theoretical predictions, the extension of the calculation methods
from $\theta=0$ to $\theta\neq 0$ is somewhat non-trivial, even in the
Hamiltonian formalism.

We performed the DMRG calculations and obtained the mass spectra mainly by the
two types of calculation methods, namely the one-point-function scheme and the
dispersion-relation scheme with several improvements to the ones originally
proposed in ref. Itou:2023img . We found that we can precisely calculate the
mass spectra of stable mesons, and the results of these two schemes are
consistent with each other even in the $\theta\neq 0$ regime. Furthermore, our
results indicate that the $\theta$-dependence of the pion mass agrees with the
analytic result, $|\cos(\theta/2)|^{2/3}$, and the sigma meson is $\sqrt{3}$
times heavier than the pion as predicted. We also confirmed the nearly
conformal behavior at $\theta=\pi$ by comparing the one-point functions of the
conformal theory on a finite interval with its analytic calculation.

This paper is organized as follows. In section 2, we revisit the continuum
$2$-flavor Schwinger model and its bosonization focusing on the low-energy
mass spectrum. We review the three methods for computing the mass spectrum in
section 3 and discuss the extension to $\theta\neq 0$ in section 4. In section
5, we present our simulation results of the mass spectra by the one-point and
dispersion-relation schemes. In section 6, we show the one-point functions at
$\theta=\pi$ and compare them with the analytic results. Section 7 is devoted
to the conclusion and discussion. In appendix A, we explain the analytic
calculation of the one-point function in the WZW model on a finite interval.
Appendix B shows the definition of the observables used in the analysis. In
appendix C, we summarize the results of the correlation-function scheme for
reference. In appendix D, we show the result of the one-point-function scheme
with a different setup from the main part. In appendix E, we discuss the
behavior of the one-point and correlation functions of the unstable eta meson.

##  2 Review of the 2-flavor Schwinger model with the $\theta$ term

In this section, we review the basic properties of the $2$-flavor Schwinger
model 222 For the original reference, see ref. Coleman:1976uz , and section 2
of our previous paper Itou:2023img  would also be useful.. The Lagrangian of
the massive $2$-flavor Schwinger model 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}^{2}\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, $\theta$ is the vacuum angle describing
the background electric flux, $m$ is the flavor-symmetric fermion mass, and
the index $f$ labels the flavor. We note that the mass dimension is
$[g]=[m]=1$, and we focus on the strongly-coupled regime, $0<m\ll g$.

###  2.1 Composite particles and global symmetry at $\theta=0$

The Coulomb interaction in $(1+1)$d is the linear potential. Therefore, the
massive multi-flavor Schwinger model is an example of confining gauge
theories, and its particle spectrum for $0<m\ll g$ is analogous to the meson
spectrum of $(3+1)$d QCD.

Let us first summarize the physics at $\theta=0$. For nonzero fermion mass
$m>0$, the $2$-flavor Schwinger model has the mass gap and the asymptotic
particle states are well defined. The global symmetry provides the quantum
numbers for their classification, and the global symmetry at $\theta=0$ is
given by

$\frac{\mathrm{SU}(2)_{V}}{\mathbb{Z}_{2}}\times(\mathbb{Z}_{2})_{G}.$ (2)

Here, $\mathrm{SU}(2)_{V}/\mathbb{Z}_{2}$ is the isospin symmetry, which acts
on the fermions as a vector-like symmetry,
$\psi\to\exp(i\varepsilon_{a}\tau_{a})\psi$ and
$\bar{\psi}\to\bar{\psi}\exp(-i\varepsilon_{a}\tau_{a})$, and we denote its
generators as $J_{a}$,

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

The $G$-parity $(\mathbb{Z}_{2})_{G}$ is defined by the combination of the
charge conjugation and the $\pi$ rotation in the isospin space,

$G=Ce^{i\pi J_{y}}.$ (4)

In addition to these internal symmetries, we also have the parity quantum
number $P$.

When $\theta=0$, there are three types of the stable particles Coleman:1976uz
. By analogy to QCD, we call them pion, sigma meson, and eta meson, and they
correspond to the lowest energy states for the following composite operators:

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

Here, $J^{PG}$ show their isospin, parity, and $G$-parity quantum numbers.
When there is no room for confusion, we abbreviate the pions $\pi_{a}$ as
$\pi$ for simplicity. For $m\ll g$, the masses of the pion and sigma meson
scale as $(m/g)^{2/3}$, and they have the specific ratio
$M_{\sigma}/M_{\pi}=\sqrt{3}$ Coleman:1976uz . Thus, the decay process
$\sigma\to\pi\pi$ is energetically forbidden, and the sigma meson is a stable
particle, unlike the case of $(3+1)$d QCD. On the other hand, the eta meson
has a heavier mass, $M_{\eta}\sim\mu=\sqrt{\frac{2}{\pi}}g$, compared with the
pion and sigma meson, but it is stable due to the parity and $G$-parity
quantum numbers. In the next subsection, we discuss how the $\theta$ term
changes the physics discussed here.

###  2.2 The particle spectrum at $\theta\neq 0$

The $\theta$ term has an important effect on the low-energy physics. The
structure of the internal global symmetry is strongly affected by the value of
$\theta$, and it can be summarized as

$\frac{\mathrm{SU}(2)_{V}}{\mathbb{Z}_{2}}\times\left\\{\begin{array}[]{ll}(\mathbb{Z}_{2})_{G}&\quad(\theta=0),\\\
\mathrm{nothing}&\quad(\theta\not=0,\pi),\\\
(\mathbb{Z}_{2})_{\mathrm{chiral}+G}&\quad(\theta=\pi).\end{array}\right.$ (8)

The $\theta$ term violates the $G$-parity explicitly as it is odd under the
charge conjugation, which makes the physics at $\theta=0$ special. We note
that $\theta=\pi$ is another special point: Although the $G$-parity itself is
violated,
$(\mathbb{Z}_{2})_{\mathrm{chiral}+G}\subset\mathrm{SU}(2)_{L}\times(\mathbb{Z}_{2})_{G}$
is kept intact, where $\mathrm{SU}(2)_{L}$ is the chiral isospin
transformation acting only on the left-moving fermions. In this section, we
consider the case $0<\theta<\pi$ and postpone our discussion for $\theta=\pi$
to the next subsection.

The Abelian bosonization gives the useful tool to analyze this system when
$m\ll g$ Coleman:1976uz , and the effective Lagrangian after integrating out
the $U(1)$ gauge field becomes

$\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],$
(9)

where $\mu^{2}=2g^{2}/\pi$ and $C=e^{\gamma}/(2\pi)$, and the symbol
$N_{\rho}[\cdot]$ denotes the normal ordering at the scale $\rho$
Coleman:1974bu . Here $\eta$ is the non-compact scalar, and $\varphi$ is the
$2\pi$-periodic scalar. In this description, the isospin
$\mathrm{SU}(2)_{V}/\mathbb{Z}_{2}$ symmetry is not manifest in the
Lagrangian, but it is realized as the symmetry acting on solitons at the
quantum level. The $G$-parity is given by $\eta\to-\eta$, and we can
explicitly see that it is a good symmetry only at $\theta=0$. At $\theta=\pi$,
the $G$-parity $\eta\to-\eta$ associated with $\varphi\to\varphi+\pi$ is the
symmetry of this Lagrangian, and this is
$(\mathbb{Z}_{2})_{\mathrm{chiral}+G}$. We note that the parity $P$ is also
violated for nonzero $\theta$, but $P^{\prime}=P\times G$ is a good quantum
number for any values of $\theta$.

Since the $\eta$ field has the mass $\mu$ at the tree level, we can integrate
it out to obtain the effective theory for the $2\pi$-periodic field $\varphi$.
In the leading order, this gives the sine-Gordon theory, and we can analyze it
by applying the optimized perturbation and the WKB analysis Coleman:1976uz .
The lightest particle state turns out to be an iso-triplet pseudoscalar, which
we identify as the pion, and its mass $M_{\pi}(\theta)$ is predicted to behave
as

$M_{\pi}(\theta)\sim\left|m\sqrt{\mu}\cos\frac{\theta}{2}\right|^{\frac{2}{3}}.$
(10)

The next lightest state is an iso-singlet scalar, and we identify it as the
$\sigma$ meson. The WKB analysis gives a specific mass relation between the
pion and sigma-meson masses,

$M_{\sigma}(\theta)=\sqrt{3}M_{\pi}(\theta).$ (11)

Up to this point, the physics is the same between $\theta=0$ and nonzero
$\theta$ except that the pion and sigma-meson masses have the nontrivial
$\theta$ dependence.

Let us move to the discussion on $\eta$. The computation of the self-energy
for $\eta$ field is not so straightforward due to the potential infrared
divergence Coleman:1976uz . We can still understand some heuristic behaviors
of the eta meson from the symmetry and the effective Lagrangian. From the
viewpoint of symmetry, the $\theta$ term violates both the parity and the
$G$-parity, and thus $\eta$ is no longer distinguished from $\sigma$ or
$\pi\pi$ scattering states by quantum numbers. This suggests the presence of
the $\eta$-$\sigma$ mixing, $\eta\leftrightarrow\sigma$, and also the opening
of the decay channel, $\eta\to\pi\pi$. Indeed, the interaction term of (9) for
$\eta$ gives

$\cos\left(\sqrt{2\pi}\eta-\frac{\theta}{2}\right)\cos\varphi=\left(\cos\frac{\theta}{2}+\sqrt{2\pi}\eta\sin\frac{\theta}{2}+\cdots\right)\cos\varphi,$
(12)

which contains the $\eta\cos\varphi$ vertex at $\theta\neq 0$. Since the
$\cos\varphi$ vertex creates the one-particle state of $\sigma$ and also the
s-wave scattering states of $\pi\pi$, this gives the microscopic explanation
of the $\eta$-$\sigma$ mixing and the $\eta\to\pi\pi$ decay.

The scalar fields $\eta$ and $\varphi$ in (9) are related to the fermion
bilinear operators as

$\bar{\psi}_{f}\frac{1\pm\gamma_{5}}{2}\psi_{f}=-\frac{1}{2}C\rho
N_{\rho}\left[\exp\left(\pm
i\left(\sqrt{2\pi}\eta-\frac{\theta}{2}-(-1)^{f}\varphi\right)\right)\right],$
(13)

for each flavor $f=1,2$. This suggests that the meson operators, (5) - (7),
should be rotated by the axial transformation to obtain the correct one-
particle states, and we find

$\displaystyle\pi_{a}$
$\displaystyle=-i\bar{\psi}\exp({i\frac{\theta}{2}\gamma^{5}})\gamma^{5}\tau_{a}\psi,$
(14) $\displaystyle\sigma$
$\displaystyle=\bar{\psi}\exp[i\left(\frac{\theta}{2}+\omega(\theta)\right)\gamma^{5}]\psi,$
(15) $\displaystyle\eta$
$\displaystyle=-i\bar{\psi}\exp[i\left(\frac{\theta}{2}+\omega(\theta)\right)\gamma^{5}]\gamma^{5}\psi.$
(16)

For the sigma and eta mesons, the finite $m$ effect gives the mixing as we
have seen in (12), and the extra diagonalization is required, which is
indicated by $\omega(\theta)$. The optimized perturbation suggests that
$\omega(\theta)$ is of the order of
$(m/g)^{4/3}\sin(\theta/2)\cos^{1/3}(\theta/2)$ at the leading order, so this
extra rotation vanishes at $\theta=0$ and $\theta=\pi$. For the pion, there is
no mixing counterpart due to the absence of the iso-triplet scalar particle,
and thus such an extra rotation would not exist. We will numerically confirm
this point in section 5.2.1.

###  2.3 Nearly conformal behaviors at $\theta=\pi$

Finally, let us consider the case of $\theta=\pi$. When we set $m=0$ in (9),
we obtain the self-dual compact boson as the low-energy effective theory,
which describes the $SU(2)_{1}$ WZW model. This is a conformal field theory
that has the central charge $c=1$ and enjoys the
$[\mathrm{SU}(2)_{L}\times\mathrm{SU}(2)_{R}]/\mathbb{Z}_{2}$ chiral symmetry.
Turning on the small $m>0$ at $\theta=\pi$ gives the symmetry-breaking
perturbation that preserves
$[\mathrm{SU}(2)_{V}/\mathbb{Z}_{2}]\times(\mathbb{Z}_{2})_{\mathrm{chiral}+G}$,
and the possible lowest-dimensional scalar operator is given by $J_{L}J_{R}$.
Here, $J_{L}$ and $J_{R}$ are holomorphic and anti-holomorphic parts of
$SU(2)\times SU(2)$ chiral symmetry generators, and $J_{L}J_{R}$ is the scalar
operator with the scaling dimension $2$. This $J_{L}J_{R}$ operator is non-
integrably marginal, and thus it can be either marginally relevant or
irrelevant depending on the sign of the coupling.333For the $2$d
$\mathbb{C}P^{1}$ sigma model at $\theta=\pi$, the sign of the coupling is on
the marginally irrelevant side and it flows back to the $SU(2)_{1}$ WZW fixed
point logarithmically according to the Haldane conjecture Haldane:1983ru ;
Affleck:1987ch . For the $2$-flavor Schwinger model at $\theta=\pi$, we can
determine its sign by going beyond the leading-order analysis when integrating
out $\eta$, and it turns out to be on the marginally relevant side
Coleman:1976uz . As a result, the system acquires the exponentially small mass
gap $\sim e^{-\\#g^{2}/m^{2}}g$, and the ground states are doubly degenerate
due to the spontaneous breaking of $(\mathbb{Z}_{2})_{\mathrm{chiral}+G}$
Dempsey:2023gib .

In the numerical analysis, we put the system in a finite interval of size $L$.
This introduces the typical energy scale of $O(\pi/L)$, and thus it is hard to
detect a tiny mass gap unless we can take an exponentially large system size
$L\gtrsim e^{\\#g^{2}/m^{2}}/g$. We expect that the model at $\theta=\pi$
looks almost like a gapless system in our simulation setup, and thus the
$SU(2)_{1}$ WZW model gives a good approximation. Here, we summarize the
analytic results of the one-point function with this approximation. The
details of the calculation are summarized in appendix A.

Let us consider the open interval, $0\leq x\leq L$, and assume that we put the
Dirichlet boundary condition in the bosonized description,

$\eta(x=0,L)=0,\quad\varphi(x=0,L)=0.$ (17)

This can be translated to the boundary condition for fermions by using (13).
We then find that

$\displaystyle\Braket{\sigma(x)}=\langle\overline{\psi}i\gamma_{5}\psi(x)\rangle=-\sqrt{\frac{e^{\gamma}\mu}{\pi
L}}\frac{1}{\sqrt{\sin(\pi x/L)}},$ (18)
$\displaystyle\Braket{\pi(x)}=\langle\bar{\psi}\tau_{3}\psi(x)\rangle=0.$ (19)

The result $\Braket{\pi(x)}=0$ follows immediately from the fact that the
boundary condition is the isospin symmetric.

We can easily generalize the result to the isospin-violating Dirichlet
boundary condition by introducing the twist angle $\Delta$. For example, if we
take

$\varphi(x=0,L)=-\Delta,$ (20)

the analytic forms of the sigma and pion one-point functions are given by

$\displaystyle\Braket{\sigma(x)}=-\sqrt{\frac{e^{\gamma}\mu}{\pi
L}}\frac{\cos\Delta}{\sqrt{\sin(\pi x/L)}},$ (21)
$\displaystyle\Braket{\pi(x)}=\sqrt{\frac{e^{\gamma}\mu}{\pi
L}}\frac{\sin\Delta}{\sqrt{\sin(\pi x/L)}}.$ (22)

As another isospin-violating Dirichlet boundary condition, we may take

$\varphi(x=0)=-\varphi(x=L)=-\Delta,$ (23)

and we here assume $|\Delta|<\pi/2$. Then, we obtain

$\displaystyle\Braket{\sigma(x)}=-\sqrt{\frac{e^{\gamma}\mu}{\pi
L}}\frac{\cos[\Delta(1-2x/L)]}{\sqrt{\sin(\pi x/L)}},$ (24)
$\displaystyle\Braket{\pi(x)}=\sqrt{\frac{e^{\gamma}\mu}{\pi
L}}\frac{\sin[\Delta(1-2x/L)]}{\sqrt{\sin(\pi x/L)}}.$ (25)

The result has a nontrivial dependence on the angle $\Delta$.

##  3 Basic calculation strategy and review at $\theta=0$

Here, we address the basics of calculation and explain our strategy to obtain
the mass spectra in the previous work at $\theta=0$ Itou:2023img . First, we
define the Hamiltonian of the $2$-flavor Schwinger model and its lattice
regularization and explain how to generate the ground state in the DMRG
method. After that, we quickly review three computational schemes for the
meson spectra, which were demonstrated in the $2$-flavor Schwinger model at
$\theta=0$ in our previous work, and summarize each pros and cons found in the
actual calculations.

###  3.1 Lattice Hamiltonian simulation of the 2-flavor Schwinger model with
the open boundary condition

We use the lattice Hamiltonian formalism to perform the numerical computation
of the $2$-flavor Schwinger model (1). As we use the same formalism with our
previous paper Itou:2023img , we keep our explanation brief when it has an
overlap.

We put the system on the open interval $0\leq x\leq L$, and perform the
canonical quantization in the temporal gauge. In this setup, we can solve the
Gauss law constraint explicitly, and the Hamiltonian can be completely
described by fermions with non-local interactions. We take the lattice
discretization of this fermionic system using the staggered fermion
Kogut:1974ag ; Susskind:1976jm  and obtain the quantum system with the finite-
dimensional Hilbert space. Let $\chi_{f,n}$ be the staggered fermion on the
lattice sites $n=0,1,\cdots,N-1$ with the lattice spacing $a$ (i.e.
$(N-1)a=L$), and then it is translated to the continuum Dirac fermion
$\psi_{f}(x)$ at $x=na$ as

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

We use this dictionary to obtain the lattice expression for various fermion
bilinear operators, and we give explicit forms for important ones in appendix
B. After these manipulations, the lattice Hamiltonian of $N_{f}$-flavor
Schwinger model consists of the gauge part $H_{J}$, the fermion kinetic term
$H_{w}$, and the mass term $H_{m}$, i.e. $H=H_{J}+H_{w}+H_{m}$, and each term
is given by

$\displaystyle H_{J}$
$\displaystyle=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},$
(27) $\displaystyle H_{w}$
$\displaystyle=-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),$
(28) $\displaystyle H_{m}$
$\displaystyle=m_{\mathrm{lat}}\sum_{f=1}^{N_{f}}\sum_{n=0}^{N-1}(-1)^{n}\chi_{f,n}^{\dagger}\chi_{f,n},$
(29)

where $J=g^{2}a/2$ and $w=1/2a$. We relate the lattice fermion mass
$m_{\mathrm{lat}}$ and the mass $m$ of the continuum theory as

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

so that we maintain the $\mathbb{Z}_{2}$ discrete chiral symmetry in the
chiral limit and prevent the additive renormalization Dempsey:2022nys .

###  3.2 Generation of the ground state

When performing the actual computation, we further map this fermionic system
to the spin system via the Jordan-Winger transformation, and its details are
explained in section 3 and appendix A of Itou:2023img  with the same notation.
After this mapping, we apply the density-matrix renormalization group (DMRG)
White:1992zz ; White:1993zza ; Schollw_ck_2005 ; Schollw_ck_2011 .

The DMRG is a variational algorithm using the matrix product state (MPS) as an
ansatz of the wave function. It can find the ground state, which minimizes the
cost function, namely the expectation value of the Hamiltonian. There is a
cutoff parameter $\varepsilon$ controlling the truncation error of the
singular-value decomposition (SVD), which determines the bond dimension $D$ of
MPS. Smaller $\varepsilon$ yields a better approximation, whereas it requires
a larger bond dimension and a higher computational cost. We choose the Néel
state as an initial state of the DMRG, which corresponds to a zero-particle
state in terms of the Schwinger model. We also impose the conservation of the
total electric charge $Q$ during the DMRG, so that the resulting MPS satisfies
$Q=0$ as required by the Gauss law.

Let us comment on the computational cost of DMRG when applying it to the
$2$-flavor Schwinger model at $\theta\neq 0$. The potential difficulty comes
from the increase of the bond dimension, because the gap of the system, namely
the pion mass, decreases as $\theta$ approaches $\pi$. At $\theta=\pi$, the
system becomes almost gapless, and the entanglement entropy is no longer
constant but increases with the lattice size $N$ as
$S_{\mathrm{EE}}\sim(c/3)\log N$ with the central charge $c=1$ in the
continuum limit. Correspondingly, the required bond dimension $D$ is expected
to increase as $D\sim N^{c/3}$. Figure 1 depicts the example plot of the bond
dimension $D$ of the MPS for the ground state as a function of $N$. Here, we
set the theoretical parameter to $m/g=0.10$ and the lattice spacing to
$a=0.25$, which are the same as the ones in the actual calculation in this
study. While $D$ is saturated to a constant value in the gapped cases of
$0\leq\theta<\pi$, it increases with $N$ at $\theta=\pi$. By fitting $D$ at
$\theta=\pi$ with $c_{1}N^{1/3}+c_{2}$, we obtain $c_{1}=288(3)$ and
$c_{2}=-563(16)$. Note that a finite spacing effect can contribute to a
discrepancy between the data and the fitting result. In this study, we take
$N\approx 300$, then the necessary bond dimension is $D\lesssim 1400$. It is a
doable calculation using a PC cluster or a single-node supercomputer.

Figure 1:  The bond dimension $D$ of the MPS is plotted against the system
size $N$ for $\theta/2\pi=0.0$, $0.3$, and $0.5$. The lattice spacing is set
to $a=0.25$. The data for $\theta/2\pi=0.5$ are fitted by
$c_{1}N^{1/3}+c_{2}$, where the dotted curve denotes the fitting result.

###  3.3 Review of computational schemes for the meson spectra at $\theta=0$

In our previous work Itou:2023img , we developed three distinct methods to
compute the mass spectrum in the Hamiltonian formalism. We demonstrated the
methods in the $2$-flavor Schwinger model at $\theta=0$ using DMRG and
confirmed that the results are consistent with each other. Let us briefly
summarize the features of the three methods.

####  3.3.1 Correlation-function scheme

The first one, the correlation-function scheme, is to employ the technique in
the Lagrangian formalism, where we obtain the meson mass from the spatial two-
point correlation function for desired composite states. In Euclidean path
integral, the meson mass can be calculated from the imaginary-time two-point
correlation function, and it is equivalent to the spatial two-point
correlation function thanks to the Lorentz invariance: When an operator $O(r)$
creates the one-particle state with the mass $M$, the asymptotic behavior of
the spatial correlator takes the Yukawa-type form,

$\displaystyle\langle\mathcal{O}(r)\mathcal{O}(0)\rangle-\langle
O(r)\rangle\langle O(0)\rangle\sim\frac{1}{r^{\alpha}}\exp(-Mr).$ (31)

Strictly speaking, the imaginary-time correlator and the spatial correlator
are no longer equivalent in the lattice regularized Hamiltonian formalism, and
the equivalence is recovered after taking the continuum limit.

While this method can apply to various models in any dimension, there is a
subtlety when using the tensor network. In eq. (31), there is the power-like
behavior $1/r^{\alpha}$ in front of the exponential factor, where
$\alpha=(d-1)/2$ for the free scalar with the spacetime dimension $d$.
However, any tensor-network-based algorithms are a kind of finite-rank
approximation of the transfer matrix, and reproducing such a power-law factor
requires a large bond dimension. Thus, we need to increase the bond dimension
until the correct asymptotic behavior is observed.

Let us point out another issue when we take the open boundary condition. The
functional form (31) assumes the translational invariance, but the presence of
boundaries violates it. Therefore, we need to take a sufficiently large system
size so that the insertion of local operators is very far from the boundaries
to make their effect negligible. When we consider the situation of the tiny
mass gap, this is difficult to achieve, and we have to take careful
extrapolations to find reliable results. In this paper, we find it difficult
to obtain the mass spectra at larger values of $\theta$, so we put the
analysis of the correlation-function scheme in appendix C instead of the main
text.

On the other hand, the correlation function still provides the most systematic
way to analyze the mixing of states. By considering the set of local operators
$O_{i}(r)$, we can compute the matrix of correlation function $\langle
O_{i}(r)O_{j}(0)\rangle$, and we may extract the information of excited
states, state-mixing, etc. by its diagonalization. This is the unique feature
of the correlation-function scheme, which is absent in the other two methods.

####  3.3.2 One-point-function scheme

The second method, the one-point-function scheme, utilizes the boundary or
defect as a source for excitations from the thermodynamic ground state. Let us
again assume that the Lorentz symmetry is approximately recovered, then we may
regard the spatial direction as the imaginary-time direction. Then, the
boundary gives some state, $|\mathrm{Bdry}\rangle$, and the one-point function
can be interpreted as

$\langle O(x)\rangle=\langle\mathrm{Bdry}|e^{-Hx}O|\Omega\rangle,$ (32)

where $|\Omega\rangle$ is the ground state and we normalize the ground-state
energy to be $0$. As $|\mathrm{Bdry}\rangle$ is translational invariant, it is
the zero mode of the momentum operator, and the asymptotic behavior of the
one-point function for the gapped system should behave as

$\displaystyle\langle\mathcal{O}(x)\rangle\sim e^{-Mx},$ (33)

where $M$ is the mass of the lightest meson with the quantum number of $O$.
This is not the Yukawa-type but the simple exponential decay unlike eq. (31),
and one may understand it with the analogy of the wall source method in
lattice QCD.

We often need to choose the symmetry-violating boundary condition when we
measure the meson mass with the nontrivial quantum number. Otherwise, the one-
point function vanishes due to the symmetry, and we cannot measure the mass of
the target particle. For example, in the previous work Itou:2023img , we used
the edge mode of the SPT state as a source of the triplet mesons by setting
$\theta=2\pi$, because the naive open boundary condition at $\theta=0$ is
isospin singlet and cannot be used to measure the pion mass.

Let us comment that the simple exponential form in (33) has nice compatibility
with the tensor network algorithms. We find that the one-point function is
insensitive to the bond dimension, which allows us to set a larger
$\varepsilon$ than that in the correlation-function scheme. However, one
cannot access the information of the off-diagonal correlators unlike the
correlation-function scheme, so this method is specialized to the lightest
particle spectrum in a given quantum number.

####  3.3.3 Dispersion-relation scheme

The third method, the dispersion-relation scheme, is the distinctive strategy
of Hamiltonian formalism. We generated the low-energy excited states as well
as the ground state. Suppose we already have the lower energy eigenstates
$\ket{\Psi_{\ell^{\prime}}}$ for the level $\ell^{\prime}=0,1,\cdots,\ell-1$.
Then the $\ell$-th excited state $\ket{\Psi_{\ell}}$ can be obtained by DMRG
with the Hamiltonian modified as

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

where the parameter $W>0$ is a weight to impose the orthogonality Wall_2012 ;
Banuls:2013jaa .

The type of the meson was identified by the quantum numbers, such as the
isospin and $G$-parity. After identifying the state of each meson by the
quantum numbers, the mass is obtained by the dispersion relation,
$E=\sqrt{M^{2}+K^{2}}$, for each meson by measuring the energy ($E$) and
momentum square ($K^{2}$). This is a basic idea of this scheme.

As we work on the open boundary, the translational invariance is explicitly
broken and the momentum operator does not provide a good quantum number. As an
ad hoc substitute, we identify the momentum square of the $\ell$-th state by
subtracting the ground-state contribution, $\Delta
K^{2}_{\ell}:=\langle\Psi_{\ell}|K^{2}|\Psi_{\ell}\rangle-\langle\Psi_{0}|K^{2}|\Psi_{0}\rangle$,
and it works well empirically.

This method is a heuristic approach and allows the detection of states other
than single mesons. This paper does not address such states, but they are
considered to be scattering states Harada:1993va . On the other hand, the
efficiency of this method depends on the way of obtaining the excited states
of interest. If we are interested in the high-level state, the computational
cost of finding the target state from the others can be a bottleneck.

##  4 Caluclation strategy for $\theta\neq 0$ regime

At $\theta\neq 0$, the ground state can be obtained by DMRG, but the expected
physics are different from the one at $\theta=0$ as discussed in sections 2.2
and 2.3. For instance, there is an operator mixing between the scalar and
pseudo-scalar in the state, and some quantum numbers are no longer exact. We
have to resolve the mixing between them to obtain the mass eigenstate.
Furthermore, a nearly conformal theory emerges at $\theta=\pi$ where the
functional form of the one-point function must be changed as shown in section
2.3. Therefore, some methodological improvements are needed.

From these situations, we take two types of calculation and introduce some
technical improvements at $\theta\neq 0$ to obtain the promising result of the
mass spectrum. The first method is an improved one-point-function scheme
combining the correlation-function and one-point-function schemes at
$\theta=0$. In this scheme, we introduce an alternative boundary condition to
control the boundary effects. The second one is the dispersion-relation
scheme, which is the same as the one at $\theta=0$. However, the parity and
$G$-parity are no longer the quantum numbers of the mesons for $\theta\neq 0$,
since the $\theta$ term explicitly breaks these symmetries. Moreover, the mass
spectra are getting lighter in the large $\theta$ regime, so that the
computational cost to obtain the excited states increases. Therefore, we also
improved the technique to generate excited states in the DMRG calculation.

###  4.1 Improved one-point-function scheme

The basic idea of this scheme is the following: First of all, we measure the
two-point correlation function for the $2\times 2$ matrix of the scalar and
pseudo-scalar operators and find the mixing angle between them by
diagonalizing the matrix. Then, we evaluate the one-point function of the bulk
operator with the mass eigenstates. Technically, to obtain the precise value
of the mass, we introduce some modified boundary conditions by twisting
fermion masses.

####  4.1.1 Determination of the mixing matrix

To resolve the operator mixing, we measure the two-point correlation function
for the $2\times 2$ matrix,

$\bm{C}_{\pm}(x,y)=\begin{pmatrix}\Braket{S_{\pm}(x)S_{\pm}(y)}_{c}&\Braket{S_{\pm}(x)PS_{\pm}(y)}_{c}\\\
\Braket{PS_{\pm}(x)S_{\pm}(y)}_{c}&\Braket{PS_{\pm}(x)PS_{\pm}(y)}_{c}\end{pmatrix},$
(35)

for the iso-triplet and iso-singlet sectors, where the subscript $c$ means the
connected correlator,
$\Braket{\mathcal{O}_{1}\mathcal{O}_{2}}_{c}=\Braket{\mathcal{O}_{1}\mathcal{O}_{2}}-\Braket{\mathcal{O}_{1}}\Braket{\mathcal{O}_{2}}$.
Here, we introduce the lattice version of the iso-triplet and iso-singlet
fermion bilinear operators: The $\tau_{3}$-components of the iso-triplet one
are denoted as

$S_{-}(x)\leftrightarrow\overline{\psi}\tau_{3}\psi(x),\quad
PS_{-}(x)\leftrightarrow-i\overline{\psi}\gamma^{5}\tau_{3}\psi(x),$ (36)

with $x=na$, and the iso-singlet ones are denoted as

$S_{+}(x)\leftrightarrow\overline{\psi}\psi(x),\quad
PS_{+}(x)\leftrightarrow-i\overline{\psi}\gamma^{5}\psi(x).$ (37)

The actual expression for the lattice fermion is obtained by using the
dictionary (26) and taking the three-lattice-point average to smear the
oscillation coming from the use of staggered fermions.

Since $S_{\pm}(x)$ and $PS_{\pm}(x)$ are Hermitian operators,
$\bm{C}_{\pm}(x,y)$ is a real symmetric matrix, which can be diagonalized by
an orthogonal matrix,

$R(\delta)=\begin{pmatrix}\cos\delta&-\sin\delta\\\
\sin\delta&\cos\delta\end{pmatrix}.$ (38)

Here the argument $\delta$ corresponds to the mixing angle, and the
eigenvalues of $\bm{C}_{\pm}(x,y)$ are the connected correlation functions of
the meson operators. Thus, we define the operators $\sigma(x)$ and $\eta(x)$
by

$\bm{C}_{+}(x,y)=R(\delta_{+})^{\mathrm{T}}\begin{pmatrix}\Braket{\sigma(x)\sigma(y)}_{c}&0\\\
0&\Braket{\eta(x)\eta(y)}_{c}\end{pmatrix}R(\delta_{+}),$ (39)

and the operator $\pi(x)$ by

$\bm{C}_{-}(x,y)=R(\delta_{-})^{\mathrm{T}}\begin{pmatrix}\hphantom{\pi(x)}*\hphantom{\pi(y)}&0\\\
0&\Braket{\pi(x)\pi(y)}_{c}\end{pmatrix}R(\delta_{-}).$ (40)

The explicit forms for $\pi(x)$, $\sigma(x)$ and $\eta(x)$ are given by

$\begin{pmatrix}*\\\
\pi(x)\end{pmatrix}=R(\delta_{-})\begin{pmatrix}S_{-}(x)\\\
PS_{-}(x)\end{pmatrix},\qquad\begin{pmatrix}\sigma(x)\\\
\eta(x)\end{pmatrix}=R(\delta_{+})\begin{pmatrix}S_{+}(x)\\\
PS_{+}(x)\end{pmatrix}.$ (41)

Note that $\delta_{\pm}$ could depend on the distance $r=|x-y|$ at which
$\bm{C}_{\pm}(x,y)$ is measured due to the boundary effect. However, we find
from actual simulation results that $\delta_{\pm}$ is almost independent of
$r$ in the bulk region, thus such a boundary effect on $\delta_{\pm}$ is
sufficiently small. We also confirm that $\delta_{\pm}$ is not sensitive to
the cutoff parameter $\varepsilon$.

Theoretically, we can estimate the mixing angle as follows. The finite-$m$
interaction term induces the presence of the $\eta$-$\sigma$ mixing. Using the
discussion around (12), we can give the analytic form for $\delta_{+}$ from
the bosonized effective Lagrangian. The $\eta$-$\sigma$ mixing is
characterized by the mass matrix,

$\mathcal{M}(\theta)=\mu^{2}\begin{pmatrix}1&A\sin(\theta/2)|\cos(\theta/2)|^{1/3}\\\
A\sin(\theta/2)|\cos(\theta/2)|^{1/3}&B|\cos(\theta/2)|^{4/3}\end{pmatrix},$
(42)

where $A$ and $B$ are $O((m/g)^{4/3})$ quantities. The real-symmetric matrix
$\mathcal{M}(\theta)$ can be diagonalized by the orthogonal matrix
$R(\omega(\theta))$ with the argument $\omega(\theta)$, which corresponds to
the extra rotation in eqs. (15) and (16). Then, $\delta_{-}\simeq\theta/2$ and
$\delta_{+}\simeq\theta/2+\omega(\theta)$ are expected by the bosonization
analysis.

####  4.1.2 Modifying the boundary and twisting fermion mass

Before measuring the bulk one-point function for each meson given in eq. (41),
let us introduce an alternative boundary condition instead of the naive open
boundary. The key idea of the one-point-function scheme is the usage of the
edge mode that appears at the boundary as the source of the meson. An
alternative boundary condition makes it possible to change the source operator
more flexibly.

Now, we consider supplemental lattice sites to both ends of the $1$d lattice.
We refer to the attached region as the wings regime 444 It comes from the term
stage wings in theater, namely the hidden parts at each side of a stage..
Then, we assign the heavier fermion mass in the wings regime than that in the
bulk, namely $m\ll m_{\rm wings}$ as shown in figure 2.

$0$$L=a(N-1)$$-aN_{\mathrm{wings}}$$L+aN_{\mathrm{wings}}$$m$$m_{\mathrm{wings}}$bulkwingswings$x$$m_{\mathrm{wings}}$
Figure 2: The schematic picture of the lattice setup with the wings regime. We
assign the mass parameter $m$ for $0\leq x\leq L$ and $m_{\mathrm{wings}}$ for
$x<0$ and $L<x$.

This setup corresponds to the Dirichlet boundary condition for the Dirac
fermion and is easily implemented by changing the mass $m$ locally in the
Hamiltonian (29). The boundary between the bulk and the wings becomes the
source of mesons. This is a basic idea of the alternative boundary condition.

In the case of the iso-singlet sector, namely the sigma and eta mesons, we
simply set the large constant, $m_{\rm wings}=m_{0}$, as the mass in the wings
regime. Similarly to the edge mode under the naive open boundary condition
that results in an iso-singlet source operator at $\theta=0$, the bulk one-
point function of the iso-singlet operators can be obtained with the simple
setup even at $\theta\neq 0$.

As for the iso-triplet sector, we apply a flavor-asymmetric twist to the
fermion mass in the wings regime 555 We utilized the property of the SPT phase
as the source of the triplet meson by shifting $\theta\rightarrow\theta+2\pi$
in our previous work. Similar analyses are shown in appendix D.. We define the
twisted mass by acting chiral rotation with the opposite phase $\pm\Delta$ for
the flavor $f=1,2$.

$0$$L=a(N-1)$$-aN_{\mathrm{wings}}$$L+aN_{\mathrm{wings}}$$m$$m_{0}e^{\pm
i\Delta\gamma^{5}}$bulkwingswings$x$$m_{0}e^{\mp i\Delta\gamma^{5}}$ Figure 3:
The schematic picture of the lattice setup with the wings regime for the iso-
triplet state calculation at $0\leq\theta<\pi$. We assign
$m_{0}e^{+i\Delta\gamma^{5}}$ ($m_{0}e^{-i\Delta\gamma^{5}}$) on the left side
of the wings and $m_{0}e^{-i\Delta\gamma^{5}}$ ($m_{0}e^{+i\Delta\gamma^{5}}$)
on the right side for the first (second) flavor.

Thus, we take

$m_{\mathrm{wings}}=\begin{cases}m_{0}\,e^{+i\Delta\gamma^{5}}&f=1,\\\
m_{0}\,e^{-i\Delta\gamma^{5}}&f=2,\end{cases}$ (43)

as the mass parameters in the wings. Since the twisted mass breaks the isospin
symmetry, the boundary between the wings and bulk can be a source of the
triplet meson. Thus, the one-point function of the pion becomes nonzero even
in the trivially gapped phase $-\pi\leq\theta\leq\pi$ by setting $\Delta\neq
0$. Note that there is a choice between having the same value of $\Delta$ on
both sides or not. Here, we take the opposite signs of $\Delta$ on both sides
of the wings as shown in figure 3 to assign the opposite directions of the
chiral rotation on each side.

Here, we address how the chiral rotation affects the lattice Hamiltonian. In
fact, the insertion of the flavor-asymmetric chiral rotation (43) also changes
the fermion kinetic term, eq. (28), of the lattice Hamiltonian. In the wings
regime, it is modified as follows:

$H_{w}=-i\sum_{f}\sum_{n\in N_{\rm
wings}}\left(w-(-1)^{n}\frac{m_{0}}{2}\sin\Delta\right)\left(\chi_{f,n}^{\dagger}\chi_{f,n+1}-\chi_{f,n+1}^{\dagger}\chi_{f,n}\right).$
(44)

As for the mass term, eq. (29), it is given by

$H_{m}=\sum_{f}\sum_{n\in N_{\rm
wings}}\left(m_{0}\cos\Delta-\frac{N_{f}g^{2}a}{8}\right)(-1)^{n}\chi_{f,n}^{\dagger}\chi_{f,n},$
(45)

taking the $O(a)$ correction (30) and the chiral rotation into account.

####  4.1.3 Improvement on the fitting function

The last step of the calculation is to fit the data of the one-point function
and extract the mass. In our previous work for $\theta=0$, we fit the data
points of $\log|\Braket{\mathcal{O}(x)}|$ by $-Mx+C$ to obtain the meson mass
$M$. However, the boundary state $|\mathrm{Bdry}\rangle$ generically contains
excited states such as two-particle states, which can give contamination. At
$\theta\neq 0$, the mass gap is expected to become smaller, and thus the
elimination of the excited-state contributions becomes more important.

To incorporate such a contribution, we change the fitting ansatz. We assume
that the one-point functions behave as

$\Braket{\mathcal{O}(x)}\sim Ae^{-Mx}+Be^{-(M+\Delta M)x},$ (46)

which is motivated by including the second-lowest state with a mass gap
$M+\Delta M$.666 Although the contribution of the two-particle states can be
Yukawa-type in general because of the integration over the relative momentum,
the exponential form turned out to be enough for the practical analysis of our
system. Then the effective mass, namely the logarithmic derivative of
$\Braket{\mathcal{O}(x)}$, is given by

$-\frac{d}{dx}\log\Braket{\mathcal{O}(x)}\sim M+\frac{\Delta M}{1+Ce^{\Delta
Mx}}.$ (47)

On the discretized lattice, we estimate the effective mass from the one-site
difference of $\Braket{\mathcal{O}(x)}$,

$-\frac{1}{a}\log\frac{\Braket{\mathcal{O}(x+a)}}{\Braket{\mathcal{O}(x)}}.$
(48)

Here, we also take the three-site average of the effective mass to suppress
the fluctuation due to the discretization with the staggered fermion.

###  4.2 Dispersion-relation scheme

The dispersion-relation scheme examines the energy eigenstates themselves and
can apply to the case of $\theta\neq 0$ straightforwardly. Even though the
parity and $G$-parity are not quantum numbers of the mesons, we can still
identify the pion and sigma meson as the lowest-energy states with the isospin
$J=1$ and $J=0$, respectively. On the other hand, the eta meson becomes
unstable at $\theta\neq 0$ as mentioned in section 2.2, and thus it should
disappear from the spectrum.

To make the dispersion relations more precise, we consider a technical
improvement in DMRG for excited states. In the DMRG with the Hamiltonian (34),
the sigma meson is relatively difficult to obtain compared with the pion. The
singlet states tend to appear at higher levels than the triplets since there
are many momentum excitations of the pion triplets at low levels. However, the
higher states generally require more computational cost and suffer from larger
systematic errors due to the lack of orthogonalities.

We deal with this problem by modifying the Hamiltonian. The point is to
generate the singlet states separately by the singlet projection, eliminating
the states with $J>0$ from the low-energy spectrum. For this purpose, we
include the isospin Casimir operator $\bm{J}^{2}$ in the Hamiltonian as the
cost term,

$H_{\ell}=H+W\sum_{\ell^{\prime}=0}^{\ell-1}\ket{\Psi_{\ell^{\prime}}}\bra{\Psi_{\ell^{\prime}}}+W_{J}\bm{J}^{2},$
(49)

where $W_{J}>0$ is a parameter to be tuned. Since $\bm{J}^{2}$ is a positive-
semidefinite operator, the last term increases the energy of the states with
nonzero isospin $J>0$. Then the condition $\bm{J}^{2}\ket{\Psi_{\ell}}=0$ is
automatically imposed on the target state $\ket{\Psi_{\ell}}$ when minimizing
the energy. Since only the singlet states are involved in DMRG, the
computational cost is reduced, which enables us to increase the number of
sweeps to improve orthogonality.

##  5 Simulation results

###  5.1 Simulation parameters

Let us explain the parameter setup in our calculation. Since the gauge
coupling $g$ has mass dimension $1$, we measure the energy scale in the unit
of $g$ by setting $g=1$. In this work, we always set the fermion mass $m=0.1$,
so the photon mass is $\mu=\sqrt{2/\pi}\simeq 0.8$. The lattice size $N$, the
lattice spacing $a$, a control parameter $\varepsilon$ of the truncation error
in SVD, and the number of iterations $N_{\mathrm{sweep}}$ in DMRG are
summarized in table 1 777 The explicit value of $a$ for the correlation-
function and one-point-function schemes is given by solving
$a\times(160-1)=0.2\times(200-1)$.. We use the C++ library of ITensor itensor
in our calculation.

scheme | $N$ | $a$ | $L=(N-1)a$ | $\varepsilon$ | $N_{\mathrm{sweep}}$
---|---|---|---|---|---
correlation function | $160$ | $0.25$ | $39.8$ | $10^{-10,-12,-14,-16}$ | $20$
one-point function | $320$ | $0.25$ | $79.9$ | $10^{-10}$ | $20$
dispersion relation | $100$ | $0.20$ | $19.8$ | $10^{-10}$ |  $50$, $80$
Table 1: Summary of the parameter setup for each calculation step. $L$ is the
physical lattice volume in the unit $g=1$.

In the calculation of the correlation function, we take a middle lattice size
$N=160$. In our previous simulation, the data of the correlation function
strongly depend on the truncation parameter $\varepsilon$, thus we take
several values of $\varepsilon$ and estimate the cutoff dependence. In the
measurement of the one-point functions, the cutoff parameter is fixed to
$\varepsilon=10^{-10}$, since we found that its cutoff dependence is
negligible. Instead, we take a large lattice size $N=320$ to see the long-
distance regime from the boundary.

To generate the ground and excited states for the dispersion-relation scheme,
we take a small lattice site $N=100$. We set $N_{\rm sweep}=50$ when we use
the cost Hamiltonian (34), where all eigenstates of the original Hamiltonian
are generated. On the other hand, we set $N_{\rm sweep}=80$ when we use the
cost Hamiltonian (49), where only singlet states are generated. Thanks to the
focus on singlet only, the number of levels to be calculated is reduced in the
latter case. Then we can increase the number of sweeps while keeping a
reasonable calculation cost 888 As we go to the higher level, the bond
dimension tends to increase. Thus, the total computational cost is not linear
on the number of the target states but higher.. The other technical simulation
parameters for each scheme will be given later.

###  5.2 Improved one-point-function scheme for $0\leq\theta<\pi$

Now, we present the results of the improved one-point-function scheme. Note
that we obtain the mass spectra of the pion and sigma meson in the range of
$0\leq\theta<\pi$. As for the eta meson, it becomes unstable in a large
$\theta$ regime. Thus, we put the analyses of the eta meson in appedix E
instead of the main text. As we show, at $\theta=\pi$, the data for the one-
point function for both pion and sigma meson do not fit the ansatz (47) since
the system becomes almost gapless and CFT-like there. We have a detailed
discussion about this point in section 6 independently.

####  5.2.1 Determination of the mixing angle

First, we resolve the mixing of the operator using the two-point correlation
function. As explained in sections 2.3 and 4.1.1, the axial rotation by
$\delta_{-}=\theta/2$ is expected for the pion while the extra rotation,
$\omega(\theta)=O((m/g)^{4/3}\sin(\theta/2)(\cos(\theta/2))^{1/3})$, is
required for the sigma and eta mesons, hence,
$\delta_{+}=\theta/2+\omega(\theta)$.

We plot the numerical results of the mixing angles $\delta_{\pm}$ as a
function of $\theta/2\pi$ in figure 4. To obtain this, we use the correlation
matrix $\bm{C}_{\pm}(x,y)$ for $x=(L-r)/2$ and $y=(L+r)/2$ with the distance
$r=15$ in our simulation.

Figure 4:  The mixing angles $\delta_{-}$ for the triplet and $\delta_{+}$ for
the singlet sector are plotted against $\theta/2\pi$. The dashed line denotes
$\theta/2$ for comparison with $\delta_{-}$. The result of fitting
$\delta_{+}$ by $\theta/2+\omega(\theta)$ is plotted as the solid curve.

The angle $\delta_{-}$ for the iso-triplet sector is shown with the empty blue
circles and mostly agrees with the line of $\theta/2$ as expected. We suspect
that the small deviation comes from the contributions of the excited state
such as the two-pion states. The angle $\delta_{+}$ for the iso-singlet sector
is shown with the empty orange squares, which have a clear deviation from
$\theta/2$. We fit the numerical result by $\theta/2+\omega(\theta)$ using eq.
(42) with parameters $A$ and $B$ and obtain the best-fit values $A=-0.23(2)$
and $B=0.76(4)$. The fitting curve is shown as the solid curve in the figure.
The success of this fitting strongly indicates the consistency with analytic
predictions on the $\eta$-$\sigma$ mixing at finite $m$.

####  5.2.2 Sigma meson masses

Now, let us show the measurement results of the one-point function for the
sigma meson and extract the mass. In the following analyses, we set the number
of lattice sites in the wings to $N_{\mathrm{wings}}=20$. The total number of
the site for the staggered fermion is then $N+2N_{\mathrm{wings}}=360$. We set
the fermion mass in the wings to $m_{\mathrm{wings}}=m_{0}=10$, which is a
hundred times larger than the mass $m=0.1$ in the bulk. We generate the ground
states of the whole system including the wings by DMRG for
$0\leq\theta/2\pi\leq 0.5$ with the interval $\Delta(\theta/2\pi)=0.1$. Since
attaching the wings does not change the bound dimension of MPS so much, the
increase in the computational cost is negligible.

The results for the one-point function $\Braket{\sigma(x)}$ of the sigma meson
operator (41) for the ground states are shown in the left panel of figure 5.

Figure 5:  (Left) The one-point function
$\log|\Braket{\sigma(x)-\sigma(L/2)}|$ of the sigma meson is plotted against
the distance $x$ from the boundary. The value at $x=L/2$ is subtracted from
$\Braket{\sigma(x)}$ to eliminate the constant shift in the bulk. The shaded
region indicates the wings. (Right) The effective mass calculated from the
one-point function in the left panel is plotted against $x$. The results of
fitting by (47) are shown as well. The fitting range is between the vertical
dashed lines.

From these data of the one-point function, we compute the effective mass and
fit the data using the improved fitting function eq. (47) in the range $5\leq
x\leq 20$ with parameters $M$, $\Delta M$, and $C$. The results are shown in
the right panel of figure 5, and the fitting results of the parameters are
summarized in table 2.

$\theta/2\pi$ | $M_{\sigma}$ | $\Delta M$ | $C$
---|---|---|---
0.0 | 0.74814(9) | 0.344(3) | 0.173(3)
0.1 | 0.7209(3) | 0.227(5) | 0.47(1)
0.2 | 0.6429(2) | 0.188(2) | 0.563(3)
0.3 | 0.5188(2) | 0.167(1) | 0.453(2)
0.4 | 0.3471(4) | 0.192(2) | 0.259(2)
Table 2:  The fitting results of the effective mass computed from the one-
point function $\Braket{\sigma(x)}$ of the sigma meson. The errors of these
values come from the fitting error.

Note that in the plot and table, we show the data except for $\theta=\pi$
since the fit ansatz, eq. (47), does not work at this $\theta$ where it closes
to the (nearly) conformal theory. We can see that the effective mass is not
completely flat but slightly curved in the bulk. It indicates that the
correction term from the excited state is not negligible 999 In our previous
work Itou:2023img , we observed small differences between the results obtained
by the three methods at $\theta=0$. A part of the differences can be explained
by neglecting the contribution from the excited state in the fitting of the
one-point function..

####  5.2.3 Pion masses

Next, let us focus on the iso-triplet meson, namely the pion. We generate the
ground state with the twisted fermion mass assigned to the wings region by
DMRG for $\theta/2\pi=0.0,0.1,\cdots,0.5$. As for the wings regime, we set the
twist angle $\Delta=0.1$, the mass $m_{0}=10$, and the size
$N_{\mathrm{wings}}=20$ in figure 3.

We measure the one-point function $\Braket{\pi(x)}$ of the pion (41), where
the rotation matrix $R_{-}$ is determined by the result in section 5.2.1. The
results of the one-point function and the corresponding effective mass are
shown in figure 6.

Figure 6:  (Left) The one-point function $\log|\Braket{\pi(x)}|$ of the pion
is plotted against the distance $x$ from the boundary. The shaded region
indicates the wings. (Right) The effective mass calculated from the one-point
function in the left panel is plotted against $x$. The results of fitting by
(47) are shown as well. The fitting range is between the vertical dashed
lines.

We obtain the nonzero one-point function thanks to the twisted mass assigned
in the wings regime. To obtain the effective mass by fitting with the ansatz
eq. (47), we utilize the data in the range $5\leq x\leq 20$ as shown in the
right panel of figure 6. The best-fit values of the parameters are summarized
in table 3.

$\theta/2\pi$ | $M_{\pi}$ | $\Delta M$ | $C$
---|---|---|---
0.0 | 0.431205(8) | 0.562(2) | 0.145(3)
0.1 | 0.418442(8) | 0.541(2) | 0.139(2)
0.2 | 0.37955(1) | 0.498(2) | 0.137(2)
0.3 | 0.31242(3) | 0.441(2) | 0.135(2)
0.4 | 0.2106(1) | 0.366(3) | 0.135(3)
Table 3:  The fitting results of the effective mass computed from the one-
point function $\Braket{\pi(x)}$ of the pion. The errors of these values come
from the fitting error.

Note that the fitting ansatz (47) does not work at $\theta=\pi$ again and we
will discuss its CFT-like behavior in section 6.

####  5.2.4 Summary of the one-point-function scheme

Let us summarize the results of the one-point-function scheme. The
$\theta$-dependent masses of the pion and sigma meson are shown in figure 7.

Figure 7:  The masses of the pion and sigma meson obtained by the one-point-
function scheme with the wings are plotted against $\theta/2\pi$. The eta-
meson mass at $\theta=0$ is also plotted for reference.

The mass of the eta meson at $\theta=0$ is also plotted for reference. The
gray dashed curve denotes $M_{\pi}(0)|\cos(\theta/2)|^{2/3}$, where the
overall coefficient $M_{\pi}(0)$ is determined by the numerical result at
$\theta=0$. The $\theta$-dependence of the pion mass is consistent with the
analytic calculation by the bosonized model in eq. (10). Indeed, we can see
that the pion mass decreases as $\theta$ increases, approaching the almost
gapless phase at $\theta=\pi$. Furthermore, the gray dash-dot curve denotes
$\sqrt{3}M_{\pi}(0)|\cos(\theta/2)|^{2/3}$. The result of the sigma-meson mass
is also consistent with the result of the WKB approximation in eq. (11). Our
result indicates that this relation holds not only around $\theta=0$ but also
in the large $\theta$ region.

###  5.3 Dispersion-relation scheme

Now, let us discuss the second calculation scheme, namely the dispersion-
relation scheme. We first generate the eigenstates up to the level $\ell=30$
using the cost Hamiltonian (34), where we set $W=10$ in our calculations. For
each state of the level $\ell$, we measure the energy gap $\Delta
E_{\ell}=E_{\ell}-E_{0}$ and the momentum square $\Delta
K_{\ell}^{2}=\Braket{K^{2}}_{\ell}-\Braket{K^{2}}_{0}$, where the result
$\Braket{K^{2}}_{0}$ of the ground state is subtracted from
$\Braket{K^{2}}_{\ell}$ to remove the nonzero offset.

The results of $\Delta E_{\ell}$ and $\Delta K_{\ell}^{2}$ for
$\theta/2\pi=0.0,0.2,0.4$ are shown in figure 8.

Figure 8:  (Top) The energy gap $\Delta E_{\ell}=E_{\ell}-E_{0}$ is plotted
against the level $\ell$ of the excited state. (Bottom) The square of total
momentum $\Delta K_{\ell}^{2}=\Braket{K^{2}}_{\ell}-\Braket{K^{2}}_{0}$ is
plotted against $\ell$ after subtracting the result of the ground state. The
left, center, and right columns correspond to $\theta/2\pi=0.0$, $0.2$, and
$0.4$, respectively. Some higher states are omitted due to insufficient
convergence.

We first identify the series of triplets with the monotonically increasing
momentum $\Delta K_{\ell}^{2}$ starting from the lowest triplet $\ell=1,2,3$
as the pions. Indeed, they have the isospin quantum numbers, $J=1$ and
$J_{z}=0,\pm 1$, which indicates that they are consistent with the quantum
numbers of the pions. On the other hand, there is another series of $J=1$
triplets at higher energy levels, which is expected to be the scattering
states investigated in ref. Harada:1993va . Furthermore, interestingly, we
also find quintet states, namely 5-fold degenerated states for
$\theta/2\pi\geq 0.3$, for instance, $\ell=23,\cdots,27$ in the right panel of
figure 8. We found that they have $\bm{J}^{2}=6$ and $J_{z}=0,\pm 1,\pm 2$.
Thus, they are the two-pion states with $J=2$, which comes from the
combination of the two pions with $J=1$.

Next, we generate only the singlet excited states up to the level $\ell=10$ by
DMRG with the Hamiltonian (49). The parameters are set to $W=10$ and $W_{J}=1$
so that DMRG could have a good convergence. Here we perform more sweeps,
$N_{\mathrm{sweep}}=80$, to improve the precision as shown in table 1.

The results of the energy gap $\Delta E_{\ell}$ and the momentum square
$\Delta K_{\ell}^{2}$ for $\theta/2\pi=0.0,0.2,0.4$ are shown in figure 9. We
first confirm that the states obtained by the Hamiltonian (49) have the
isospin quantum numbers consistent with $\bm{J}^{2}=0$ and $J_{z}=0$.
Furthermore, $\Delta E_{\ell}$ and $\Delta K_{\ell}^{2}$ in figure 9 are
almost consistent with those in figure 8, while some data show small
discrepancy. The discrepancy appears in the states, which have almost
degenerate energies with other states in figure 8. We expect that the results
obtained by the singlet projection are more reliable since they can preserve
the orthogonal condition without accidental degeneracy. Therefore, to obtain
the dispersion relation of the sigma and eta mesons, we use the results of the
singlet projection in the following.

Figure 9:  These results are obtained by DMRG with the singlet projection
(49). (Top) The energy gap $\Delta E_{\ell}=E_{\ell}-E_{0}$ is plotted against
the level $\ell$ of the excited state. (Bottom) The square of total momentum
$\Delta K_{\ell}^{2}=\Braket{K^{2}}_{\ell}-\Braket{K^{2}}_{0}$ is plotted
against $\ell$ after subtracting the result of the ground state. The left,
center, and right columns correspond to $\theta/2\pi=0.0$, $0.2$, and $0.4$,
respectively.

Now, let us focus on the results of the singlet projection in figure 9. We
first identify the lowest series of singlet states with the monotonically
increasing momentum as the sigma meson. At $\theta=0$, we can confirm that the
corresponding state has the positive $G$-parity, namely $\Braket{G}>0$. We
also find the singlet states with $\Braket{G}<0$, namely the eta meson, at
$\theta=0$. For $\theta\neq 0$, the eta meson is no longer a stable particle,
and thus the corresponding states might be replaced by some scattering state.
We find a series of eta-meson-like states at $\theta/2\pi=0.1$ but not for
$\theta/2\pi\geq 0.2$. By gradually increasing $\theta$ in the range
$0.1\leq\theta/2\pi\leq 0.2$, we observe that a state with $\Braket{G}<0$ in
the spectrum suddenly changes into a state with $\Braket{G}>0$.

Finally, we investigate the dispersion relations using the results of the pion
and sigma meson. We plot the energy gap $\Delta E_{\ell}$ against the momentum
square $\Delta K_{\ell}^{2}$ for each meson in figure 10.

Figure 10:  The energy gap $\Delta E_{\ell}$ is plotted against the square of
total momentum $\Delta K_{\ell}^{2}$ for $\theta/2\pi=0.0,0.1,\cdots,0.5$. The
states identified as the same meson are denoted by the same symbol. We fit the
data points for each meson by $\Delta E=\sqrt{b^{2}\Delta K^{2}+M^{2}}$, and
the results are shown by the broken curves. The markers with error bars at the
left endpoints denote the extrapolated value and its error from the fitting.

The states identified as the same meson are plotted by the same symbol. Then
we fit these 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$. The parameter $b$ is included taking account of possible
lattice artifacts on the momentum. As for the pion at $\theta=\pi$, we fit the
data points by $\Delta E=\sqrt{b^{2}\Delta K^{2}}$ assuming $M=0$ because the
fitting including the parameter $M$ is unstable due to the square root of a
small number.101010The error of $M$ is roughly $O(1/M)$ considering the error
propagation.

The $\theta$-dependent masses of the pion, sigma, and eta meson are summarized
in figure 11. Now, we obtain the masses for the pion and sigma meson at
$\theta=\pi$, which are almost zero as expected.

Figure 11:  The masses of the pion, sigma, and eta meson obtained by the
dispersion-relation scheme are plotted against $\theta/2\pi$.

Again, the gray dashed curve denotes the calculation in the bosonized model,
$M_{\pi}(0)|\cos(\theta/2)|^{2/3}$, where the overall coefficient $M_{\pi}(0)$
is determined by the numerical result at $\theta=0$. The gray dash-dot curve
denotes $\sqrt{3}M_{\pi}(0)|\cos(\theta/2)|^{2/3}$. Both the pion and sigma-
meson masses almost agree with the analytic predictions.

Let us make some comments on the fate of the eta meson. We find some states
similar to the eta meson, but the energies of such states become much higher
for $\theta/2\pi\geq 0.2$ than that for $\theta/2\pi\leq 0.1$. In fact, the
data points of $\Delta E_{\ell}$ and $\Delta K_{\ell}^{2}$ of these states are
no longer on a smooth line of the dispersion relation. Thus, it looks like the
eta meson has disappeared from the spectrum, which is consistent with that the
eta meson is no longer stable and decays into the pions or mixes with the
sigma meson. It also causes the relatively large error of the sigma-meson mass
at $\theta/2\pi=0.2$. Indeed, the spectrum is changed drastically around
$\theta/2\pi=0.2$. It suggests that the sigma meson is strongly contaminated
by mixing with the eta meson and they might be in the process of being
exchanged and mixed.

##  6 One-point function at $\theta=\pi$

The massive $2$-flavor Schwinger model at $\theta=\pi$ has the tiny mass gap
$M_{\pi}\sim e^{-\\#g^{2}/m^{2}}g$ Dempsey:2023gib , and we expect almost
conformal behavior. We here compare the behaviors of the one-point function
with its analytic computation on a finite interval with the $SU(2)_{1}$ WZW
conformal theory as shown in section 2.3.

For the sigma meson operator, the analytic result of the one-point function
$\Braket{\sigma(x)}$ is given by eq. (18). The corresponding numerical result
obtained by DMRG at $\theta=\pi$ is shown in the left panel of figure 12. In
this calculation, we use the lattice equipped with the wings regime, in which
the large fermion mass $m_{\mathrm{wings}}=m_{0}=10$ is assigned as in section
5.2.2.

We fit the data points of $\Braket{\sigma(x)}$ by $A/\sqrt{\sin(\pi x/L)}+B$
with the parameters $A$ and $B$. Then we obtain the best-fit values
$A=-0.079922(8)$ and $B=-0.014741(9)$. The corresponding fitting curve is
plotted in the left panel of figure 12. We find that the numerical result
agrees with (18) in the bulk region. The deviation around the boundary can be
interpreted as the contribution of the eta meson and the artifact of the
staggered fermion.111111 The local operators of the Dirac fermion consist of
the multi-point operators of the staggered fermion such as
$\chi_{n}^{\dagger}\chi_{n+1}$. Thus, the boundary condition for the Dirac
fermion is not straightforwardly translated into that for the staggered
fermion. In particular, it violates the discrete chiral symmetry explicitly as
it is realized as the one-unit lattice translation. We would like to note that
the nonzero value of the constant part $B$ is quite suggestive as it can be
regarded as the expectation value for the chiral condensate detecting the
$(\mathbb{Z}_{2})_{G+\mathrm{chiral}}$ symmetry breaking, while we need more
systematic studies to make it conclusive.

Figure 12:  The one-point functions at $\theta=\pi$ of the sigma meson
$\Braket{\sigma(x)}$ (left) and the pion $\Braket{\pi(x)}$ (right) are plotted
against the distance $x$ from the boundary. We use the lattice equipped with
the wings regime, setting $N=320$ in the bulk and $N_{\mathrm{wings}}=20$ in
the wings. The fermion mass in the wings is $m_{\mathrm{wings}}=m_{0}$ for the
sigma meson and $m_{\mathrm{wings}}=m_{0}e^{\pm i\Delta\gamma^{5}}$ for the
pion with $m_{0}=10$. In the latter case, the twist parameters are
$\Delta=+0.1$ and $-0.1$ on the left and right sides of the wings,
respectively. The fitting curves of $A/\sqrt{\sin(\pi x/L)}+B$ for the sigma
meson and $A\sin\left[\Delta(1-2x/L)\right]/\sqrt{\sin(\pi x/L)}$ for the pion
are shown as well, where the fitting range is between the vertical dashed
lines.

For the pion, we consider the lattice with the wings and take the large
fermion mass with the flavor-asymmetric chiral rotation
$m_{\mathrm{wings}}=m_{0}e^{\pm i\Delta\gamma^{5}}$ in the wings regime as
shown in figure 3. Here, we take $m_{0}=10$ and $\Delta=0.1$. The directions
of the twist are opposite, so that we obtain eq. (25) which has a nontrivial
$\Delta$ dependence.

The numerical result of the pion one-point functions $\Braket{\pi(x)}$ is
displayed in the right panel of figure 12. We fit the data points by
$A\sin\left[\Delta(1-2x/L)\right]/\sqrt{\sin(\pi x/L)}$ and obtain
$A=0.0663(1)$. As can be seen from the figure, the numerical results mostly
agree with the analytic result, eq. (25) in the bulk region. To confirm the
nontrivial $\Delta$ dependence, we also check that the one-point function
behaves as eq. (25) for other choices of the twist parameter,
$\Delta=0.2,0.3$. Therefore, it is numerically confirmed that the pion and
sigma meson in the $2$-flavor Schwinger model at $\theta=\pi$ are well
described by SU(2) level-1 WZW CFT on a finite interval.

##  7 Conclusion and discussion

Extending our previous study Itou:2023img , we investigate the physics of the
massive $2$-flavor Schwinger model at nonzero $\theta$ by employing the DMRG
to its lattice Hamiltonian formalism. In ref. Itou:2023img , we have proposed
three methods to obtain the meson spectra at $\theta=0$. At nonzero values of
$\theta$, the mass gap of the system is expected to become smaller and some
quantum numbers that distinguish mesons are no longer exact, so that we give
their several improvements in section 4 to find the precise results in such
cases.

Now, we summarize the $\theta$-dependent masses of the pion and sigma meson in
figure 13, which are obtained by the improved one-point-function and
dispersion-relation schemes. The eta-meson mass at $\theta=0$ is also plotted
for reference. As we can see in the plot, the results of the two schemes are
consistent with each other. We also compare the numerical results with the
analytic calculation in the bosonized model shown in section 2.2. The analytic
result of the pion mass (10) is depicted by the gray dashed curve in figure
13, normalized by the numerical data at $\theta=0$. We find that the numerical
result agrees with the analytic prediction. Furthermore, the sigma-meson mass
satisfies the relation, $M_{\sigma}/M_{\pi}=\sqrt{3}$, even for large $\theta$
as shown by the gray dash-dot curve. This result indicates that the
semiclassical analysis by the WKB-type approximation gives almost the correct
answer for the wide range of $\theta$. It is still interesting to clarify the
mechanism of this agreement and the range of applicability.

Figure 13:  The $\theta$-dependent meson masses obtained by the one-point-
function and dispersion-relation schemes are compared. The gray dashed curve
denotes the analytic calculation in the bosonized model,
$M_{\pi}(\theta)=M_{\pi}(0)|\cos(\theta/2)|^{2/3}$, where the overall
coefficient $M_{\pi}(0)$ is determined by averaging the results of the two
schemes at $\theta=0$. The gray dash-dot curve denotes the predicted sigma-
meson mass, $\sqrt{3}M_{\pi}(0)|\cos(\theta/2)|^{2/3}$.

At $\theta=\pi$, it is expected that the $2$-flavor Schwinger model becomes
almost CFT-like due to the exponentially small mass gap. In this case, it is
not feasible to obtain the meson masses from the one-point functions unless we
can manage the exponentially large system size. Instead, we compare the one-
point functions themselves with the calculation in the WZW model on a finite
interval. As shown in section 6, the numerical results agree with the analytic
calculation in the bulk, which provides numerical evidence that the $2$-flavor
Schwinger model at $\theta=\pi$ is well approximated by the $SU(2)_{1}$ WZW
model. In addition to the WZW behavior, we find the additional constant
contribution to the one-point function. This may give the hint of the
spontaneous breaking of the $(\mathbb{Z}_{2})_{G+\mathrm{chiral}}$ symmetry,
which requires more systematic comparisons of various volumes and also the
continuum limit to be conclusive.

The eta meson is no longer the stable particle for $\theta\neq 0$ since the
$G$-parity is broken. Indeed, the behaviors of the one-point and correlation
functions of the eta meson are quantitatively different from the other mesons.
We consider it consistent with the theoretical expectations for the
$\eta\rightarrow\pi\pi$ decay and the $\eta$-$\sigma$ mixing.

Overall, our results are consistent with the analytic predictions for the
massive $2$-flavor Schwinger model at small mass with nonzero $\theta$. We
would like to emphasize that our calculation covers the region with the large
$\theta$, which could not be accessed by the reweighting technique in the
conventional Monte Carlo method Fukaya:2003ph . This becomes possible thanks
to the efficiency of the DMRG for the lattice Hamiltonian formalism. It would
be interesting and important to extend the applicability of the Hamiltonian-
based approach to other systems, that suffer from the sign problem in the
Monte Carlo methods, such as the finite-density systems in higher dimensions.

###### Acknowledgements.

We would like to thank S. Aoki, M. Honda, K. Murakami, S. Takeda, and A. Ueda
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 E. I. was partially supported by the JSPS (Grants No. JP21H05190,
JP23H05439), and JST PRESTO Grant Number JPMJPR2113. This work was partially
supported by Japan Society for the Promotion of Science (JSPS) KAKENHI Grant
No. 23K22489 (Y.T.), and by Center for Gravitational Physics and Quantum
Information (CGPQI) at Yukawa Institute for Theoretical Physics.

##  Appendix A Compact boson on the finite interval

In this appendix, we analytically study the behaviors of various one-point
functions of the $2$-flavor Schwinger model at $\theta=\pi$ on the open
interval. The low-energy properties of the $2$-flavor Schwinger model at
$\theta=\pi$ are well described by the $\mathrm{SU}(2)$ level-$1$ Wess-Zumino-
Witten (WZW) model with the marginally relevant $J\bar{J}$ deformation. When
the size of the interval is not so large, the $J\bar{J}$ deformation is
negligible and we should be able to use the $\mathrm{SU}(2)_{1}$ WZW CFT as an
effective description. This is equivalent to the compact boson at the self-
dual radius, which is a free theory, so we can compute various correlation
functions analytically.

Within the approximations discussed above, the effective action becomes

$\mathcal{L}_{\mathrm{eff}}=\frac{1}{4\pi}(\partial_{\mu}\varphi)^{2},$ (50)

with $\varphi(x)\sim\varphi(x)+2\pi$. The mapping of operators between the
Schwinger model and the effective theory is given by

$\displaystyle\overline{\psi}\psi\leftrightarrow-\cos\frac{\theta}{2}\cos\varphi,\quad-\overline{\psi}i\gamma_{5}\psi\leftrightarrow\sin\frac{\theta}{2}\cos\varphi,$
(51)
$\displaystyle\overline{\psi}\tau_{3}\psi\leftrightarrow-\sin\frac{\theta}{2}\sin\varphi,\quad-\overline{\psi}i\gamma_{5}\tau_{3}\psi\leftrightarrow-\cos\frac{\theta}{2}\sin\varphi.$
(52)

For the open interval, we need to specify the boundary condition, and we here
use the Dirichlet boundary condition,
$\varphi(x_{1}=0,x_{2})=\varphi(x_{1}=L,x_{2})=0$. Let us then compute

$\langle\mathrm{e}^{\pm
i\varphi(x_{*})}\rangle=\frac{1}{Z}\int\mathcal{D}\varphi\exp(-\int\frac{1}{4\pi}(\partial\varphi)^{2}d^{2}x\pm
i\varphi(x_{*})).$ (53)

As we have the symmetry $\varphi\to-\varphi$ including the boundary condition,
we have
$\langle\mathrm{e}^{i\varphi(x_{*})}\rangle=\langle\mathrm{e}^{-i\varphi(x_{*})}\rangle$
so let us focus on the $+$ sign. The classical equation of motion is

$\frac{1}{2\pi}\partial^{2}\varphi(x)+i\delta(x-x_{*})=0,$ (54)

and we denote its solution as $\varphi_{\mathrm{cl}}$. Shifting $\varphi$ in
(53) as $\varphi_{\mathrm{cl}}+\varphi$, we obtain

$\displaystyle\langle e^{i\varphi(x_{*})}\rangle$
$\displaystyle=\frac{1}{Z}\int\mathcal{D}\varphi\exp\left(-\int\frac{1}{4\pi}(\partial\varphi+\partial\varphi_{\mathrm{cl}})^{2}d^{2}x+i(\varphi(x_{*})+\varphi_{\mathrm{cl}}(x_{*}))\right)$
$\displaystyle=\mathrm{e}^{-\int\frac{1}{4\pi}(\partial\varphi_{\mathrm{cl}})^{2}+i\varphi_{\mathrm{cl}}(x_{*})}$
$\displaystyle=\mathrm{e}^{i\varphi_{\mathrm{cl}}(x_{*})/2}.$ (55)

$0$$-L$$L$$2L$$\ominus$$\oplus$$\ominus$$\oplus$$\oplus$$(x_{*},0)$$(-x_{*},0)$$(2L-x_{*},0)$$(2L+x_{*},0)$$(-2L+x_{*},0)$
Figure 14: The schematic picture of the setup for the mirror-image method.

Thus, the problem reduces to find the classical solution
$\varphi_{\mathrm{cl}}$ satisfying the boundary condition, and we can use the
mirror-image method by extending the domain of $\varphi_{\mathrm{cl}}$ (see
figure 14):

$\displaystyle\frac{1}{2\pi}\partial^{2}\varphi_{\mathrm{cl}}$
$\displaystyle=-i\sum_{n=-\infty}^{\infty}\left[\delta(x-(x_{*}+2nL))-\delta(x-(-\bar{x}_{*}+2nL))\right]$
$\displaystyle=-\frac{i}{2\pi}\partial^{2}\sum_{n}\log\frac{|x-x_{*}+2nL|}{|x+\bar{x}_{*}+2nL|},$
(56)

where we use the complex coordinate. Thus,

$\varphi_{\mathrm{cl}}(x)=-i\sum_{n}\log\frac{|x-x_{*}+2nL|}{|x+\bar{x}_{*}+2nL|}=-i\log\left|\frac{\sin\left(\frac{\pi}{2L}(x-x_{*})\right)}{\sin\left(\frac{\pi}{2L}(x+\bar{x}_{*})\right)}\right|.$
(57)

Substituting this result formally, we find

$\langle
e^{i\varphi(x_{*})}\rangle=\mathrm{e}^{i\varphi_{\mathrm{cl}}(x_{*})/2}=\left(\frac{\sin\frac{\pi\epsilon}{2L}}{\sin\left(\frac{\pi}{2L}(x_{*}+\bar{x}_{*})\right)}\right)^{1/2},$
(58)

where we introduce the parameter $\epsilon$ for the point-splitting
regularization. With the normal ordering procedure at the scale $\rho$, we end
up with multiplying the factor
$\exp(\frac{1}{2}K_{0}(\rho\epsilon))\simeq\sqrt{\frac{2}{e^{\gamma}\rho\epsilon}}$
as $\epsilon\to 0$, and we get

$\langle N_{\rho}[e^{i\varphi(x_{*})}]\rangle=\sqrt{\frac{\pi}{e^{\gamma}\rho
L}}\frac{1}{\sqrt{\sin\left(\frac{\pi}{2L}(x_{*}+\bar{x}_{*})\right)}}.$ (59)

Setting $x_{*}$ to be real (i.e. the location of the source is set at
$x_{2}=0$), we have

$\langle N_{\rho}[\cos\varphi(x_{*})]\rangle=\sqrt{\frac{\pi}{e^{\gamma}\rho
L}}\frac{1}{\sqrt{\sin(\pi x_{*}/L)}},\quad\langle
N_{\rho}[\sin\varphi(x_{*})]\rangle=0.$ (60)

Therefore, for example, the expectation value of
$\overline{\psi}i\gamma_{5}\psi$ behaves as follows with the Dirichlet
boundary condition, when assuming the $SU(2)_{1}$ WZW CFT as the effective
description at $\theta=\pi$:

$\displaystyle\langle\overline{\psi}i\gamma_{5}\psi(x_{*})\rangle$
$\displaystyle=-2C\sqrt{\mu\rho}\langle N_{\rho}[\cos(\varphi(x_{*}))]\rangle$
$\displaystyle=-\sqrt{\frac{e^{\gamma}\mu}{\pi L}}\frac{1}{\sqrt{\sin(\pi
x_{*}/L)}}.$ (61)

On the other hand, we get $\langle\overline{\psi}\tau_{3}\psi\rangle=0$ in
this flavor-symmetric Dirichlet boundary condition.

##  Appendix B Definition of the observables

In this appendix, we show the explicit forms of the observables used to obtain
the mass spectrum. Let us start with the lattice version of the scalar and
pseudo-scalar operators, $S_{f,n}$ and $PS_{f,n}$, for the flavor $f$ at the
site $n$. Their explicit forms are given by rewriting the local operators
$\bar{\psi}\psi$ and $-i\bar{\psi}\gamma^{5}\psi$ in terms of the staggered
fermions (26). If we choose the same basis of the $(1+1)$d gamma matrices as
in section 3.1, we obtain

$S_{f,n}:=\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}),$
(62)
$PS_{f,n}:=\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}),$
(63)

for the lattice sites $n=1,2,\cdots,N-2$. Here we take the average of the
three neighboring sites to remove the oscillation from the staggered phase.

Next, we define some global observables to measure the isospin quantum numbers
$J_{z}$, $\bm{J}^{2}$, and the total momentum $K$. These operators are used in
the dispersion-relation scheme. Unlike the local operators defined above,
these global operators act on the whole lattice. Let us consider the isospin
operators (3). Using the staggered fermion, the lattice version of $J_{z}$ in
terms of the staggered fermions is given 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).$
(64)

It is convenient to define the lattice version of $J_{\pm}=J_{x}\pm iJ_{y}$ by

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

Then the isospin Casimir operator $\bm{J}^{2}$ can be defined as the
combination of the operators above by

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

We can evaluate the expectation values of these operators systematically using
the so-called matrix product operator.

Finally, we define the total momentum operator, which can be used to identify
the momentum excitation Banuls:2013jaa . The continuum description of the
gauge invariant momentum operator is given by

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

which commutes with the Hamiltonian of the continuum theory under the periodic
boundary condition using the Gauss-law constraint. In our case with the open
boundary condition, $K$ does not commute with the Hamiltonian, so it is no
longer the exact quantum number. However, the momentum operator is still
useful as an approximation to obtain the dispersion relation. After the gauge
fixing $U_{n}=1$, the lattice version of the momentum operator is defined by

$K=\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}).$
(69)

This operator does not exactly commute with the terms $H_{w}$ (28) and $H_{J}$
(27) of the lattice Hamiltonian due to the open boundary and the finite
lattice spacing effect.

##  Appendix C Results of the correlation-function scheme

In the Euclidean lattice gauge theory, the mass spectrum is obtained from the
correlation function in the imaginary time direction. In our previous work
focusing on $\theta=0$, we consider an analogous approach in the Hamiltonian
formalism, examining the equal-time spatial correlation function. In this
appendix, we extend this approach, the correlation-function scheme, to the
case of $\theta\neq 0$ and show the results for the stable mesons, namely the
pion and sigma meson. The correlation function of the unstable eta meson is
discussed in appendix E.2. As we mentioned in section 4.1.1, the two-point
connected correlation functions,

$C_{\pi}(r)=\Braket{\pi(x)\pi(y)}_{c},$ (70)
$C_{\sigma}(r)=\Braket{\sigma(x)\sigma(y)}_{c},$ (71)

with the distance $r=|x-y|$ are obtained as the eigenvalues of the correlation
matrix (35). To reproduce the correct asymptotic behavior of the correlation
functions, the cutoff parameter $\varepsilon$ has to be small enough in DMRG,
resulting in a large bond dimension of the MPS.

We generated the MPS of the ground state for $0\leq\theta\leq\pi$ with
$\varepsilon=10^{-10}$, $10^{-12}$, $10^{-14}$, and $10^{-16}$. The number of
lattice sites is set to $N=160$ in this scheme. Then the correlation functions
of the pion and sigma meson are computed by changing $x$ and $y$ symmetrically
as $x=(L-r)/2$ and $y=(L+r)/2$ for $0\leq r\leq L/2$. The results with
$\varepsilon=10^{-16}$ are shown in figure 15 in log scale. The change of the
gradient reflects the change of the mass depending on $\theta$.

Figure 15:  The correlation function of the pion
$\log|\Braket{\pi(x)\pi(y)}_{c}|$ (left) and the sigma meson
$\log|\Braket{\sigma(x)\sigma(y)}_{c}|$ (right) are plotted against the
distance $r=|x-y|$. The cutoff parameter in DMRG is $\varepsilon=10^{-16}$.

To determine the mass of each meson, we use the so-called effective mass,
namely the logarithmic derivative of the correlation function
$C_{\mathcal{O}}(r)$,

$M_{\mathcal{O},\mathrm{eff}}(r)=-\frac{d}{dr}\log C_{\mathcal{O}}(r).$ (72)

On the lattice, we define the effective mass by the three-point average,

$M_{\mathcal{O},\mathrm{eff}}(r):=\frac{1}{4}\tilde{M}_{\mathcal{O},\mathrm{eff}}(r-2a)+\frac{1}{2}\tilde{M}_{\mathcal{O},\mathrm{eff}}(r)+\frac{1}{4}\tilde{M}_{\mathcal{O},\mathrm{eff}}(r+2a),$
(73)

where $\tilde{M}_{\mathcal{O},\mathrm{eff}}(r)$ is a discretized logarithmic
derivative of the correlation function,

$\tilde{M}_{\mathcal{O},\mathrm{eff}}(r)=-\frac{1}{2a}\log\frac{C_{\mathcal{O}}(r+2a)}{C_{\mathcal{O}}(r)}.$
(74)

Here the factor $2a$ comes from the step size of changing $r$.

Since we consider the spatial correlation function, the leading asymptotic
behavior for $r\to\infty$ is not purely the exponential decay, but the Yukawa-
type form,

$C_{\mathcal{O}}(r)\sim\frac{1}{r^{\alpha}}\exp(-Mr).$ (75)

For example, we have $\alpha=1/2$ in the case of the $(1+1)$d free massive
boson. The corresponding effective mass for the Yukawa-type correlator is
given by

$M_{\mathcal{O},\mathrm{eff}}(r)\sim\frac{\alpha}{r}+M,$ (76)

where an $O(1/r)$ contribution exists on top of the actual meson mass $M$.
With this in mind, let us see our numerical results. The effective masses
computed from the correlation function of the pion and sigma meson are plotted
against $1/r$ in the top and bottom rows of figure 16, respectively. As
$\varepsilon$ is decreased, namely the bond dimension is increased, the
effective mass approaches $\propto 1/r$ asymptotically. This behavior is
observed for $\theta\neq 0$ as well as for $\theta=0$, although the results
are affected by the boundary when $r$ approaches the system size $L$.

Figure 16:  The effective masses $M_{\pi,\mathrm{eff}}(r)$ of the pion (top)
and $M_{\sigma,\mathrm{eff}}(r)$ of the sigma meson (bottom) are plotted
against $1/r$ for various values of $\varepsilon$. The left, center, and right
columns are the results at $\theta/2\pi=0.0$, $0.2$, and $0.4$, respectively.
The data points for $\varepsilon=10^{-16}$ are fitted by $\alpha/r+M$. The
fitting range is between the vertical dashed lines. The fitting result is
depicted by the shaded band with the systematic error.

Assuming the asymptotic behavior of Yukawa type (76), we estimate the meson
mass by the linear extrapolation $1/r\rightarrow 0$ of the effective mass
(73), which is performed by fitting the data points by $\alpha/r+M$. Here the
MPS generated with $\varepsilon=10^{-16}$ is used to suppress the effect of
truncating the bound dimension as much as possible. To evaluate the systematic
errors from the uncertainty of the fitting range, we try fitting many times by
changing the fitting range inside a given maximum region, and we obtain a
histogram of the fitting results. The best-fitting result and its error are
estimated from the position and width of the peak. We set the maximum fitting
region $0.075\leq 1/r\leq 0.15$. The fitting results for the pion and sigma
meson are summarized in table 4 and 5, respectively. The fitting curves are
depicted by the purple shadows in figure 16.

$\theta/2\pi$ | $M_{\pi}$ | $\alpha$
---|---|---
0.0 | 0.431(1) | 0.476(9)
0.1 | 0.416(1) | 0.49(1)
0.2 | 0.376(3) | 0.51(2)
0.3 | 0.306(4) | 0.55(4)
0.4 | 0.201(4) | 0.62(4)
0.5 | 0.031(4) | 0.88(4)
Table 4:  The fitting results of the effective mass $M_{\pi,\mathrm{eff}}(r)$ of the pion by $\alpha/r+M$. The errors of these values come from the systematic error from the uncertainty of the fitting range. $\theta/2\pi$ | $M_{\sigma}$ | $\alpha$
---|---|---
0.0 | 0.722(6) | 0.83(5)
0.1 | 0.699(4) | 0.79(4)
0.2 | 0.624(4) | 0.78(4)
0.3 | 0.503(3) | 0.81(3)
0.4 | 0.334(5) | 0.85(4)
0.5 | 0.134(6) | 0.64(5)
Table 5:  The fitting results of the effective mass
$M_{\sigma,\mathrm{eff}}(r)$ of the sigma meson by $\alpha/r+M$. The errors of
these values come from the systematic error from the uncertainty of the
fitting range.

The fitting of the effective mass at $\theta=\pi$ results in $M_{\pi}\sim 0$,
which indicates the correlation function of the pion is almost power-law
$C_{\pi}(r)\sim 1/r^{\alpha}$. This result is reasonable since the system
becomes CFT-like there. On the other hand, we have a relatively large mass of
the sigma meson at $\theta=\pi$ compared with the pion due to the boundary
effect (finite size effect). Indeed, the region where the $1/r$ behavior can
be observed is less clear for the sigma meson. The $\theta$-dependent meson
masses are summarized in figure 17. The result is almost consistent with that
of the one-point-function scheme in figure 7 and the dispersion-relation
scheme in figure 11. The $\theta$-dependence of the pion mass agrees with the
calculation in the bosonized model in (10), and the sigma-meson mass is close
to the result of the WKB approximation in (11).

Figure 17:  The masses of the pion and sigma meson obtained by the
correlation-function scheme are plotted against $\theta/2\pi$. The eta-meson
mass at $\theta=0$ is also plotted for reference. The gray dashed curve
denotes the analytic calculation by the bosonized model,
$M_{\pi}(0)|\cos(\theta/2)|^{2/3}$, where the overall coefficient $M_{\pi}(0)$
is determined by the numerical result at $\theta=0$. The gray dash-dot curve
denotes $\sqrt{3}M_{\pi}(0)|\cos(\theta/2)|^{2/3}$.

##  Appendix D One-point-function scheme with the original setup

In section 5.2, we show the results of the improved one-point-function scheme
on the lattice equipped with the wings. In this appendix, we show the results
of the one-point-function scheme with the naive open boundary condition used
in our previous work Itou:2023img . Then we check the consistency of these
results. Here the way of analysis is the same as section 5.2, and the
difference is only the boundary condition.

The naive open boundary can be a source of the singlet meson. We measure the
one-point function $\Braket{\sigma(x)}$ of the sigma meson (41) for the ground
state, where we use the result in section 5.2.1 to determine $R_{+}$. The
result of $\Braket{\sigma(x)}$ is shown in the left panel of figure 18, which
decays almost exponentially in the bulk region. As we did in section 5.2.2, we
compute the effective mass (48) and fit it by (47) in the range $5\leq x\leq
20$. The result is shown in the right panel of figure 18, and the fitting
parameters are summarized in table 6.

Figure 18:  (Left) The one-point function $\log|\Braket{\sigma(x)-\sigma(L/2)}|$ of the sigma meson is plotted against the distance $x$ from the boundary for $\varepsilon=10^{-10}$. The value at $x=L/2$ is subtracted from $\Braket{\sigma(x)}$ to eliminate the constant shift in the bulk. (Right) The effective mass calculated from the one-point function in the left panel is plotted against $x$. The results of fitting by (47) are shown as well. The fitting range is between the vertical dashed lines. $\theta/2\pi$ | $M_{\sigma}$ | $\Delta M$ | $C$
---|---|---|---
0.0 | 0.74840(9) | 0.317(2) | 1.194(7)
0.1 | 0.7205(1) | 0.186(1) | 1.645(8)
0.2 | 0.6414(3) | 0.154(2) | 1.45(2)
0.3 | 0.5156(8) | 0.142(3) | 0.94(4)
0.4 | 0.3453(8) | 0.178(3) | 0.65(2)
Table 6:  The fitting results of the effective mass computed from the one-
point function $\Braket{\sigma(x)}$ of the sigma meson. The errors of these
values come from the fitting error.

In section 5.2, we applied the flavor-asymmetric chiral rotation to obtain the
nonzero one-point function of the pion. Here, we instead use the edge mode of
the SPT state induced by shifting $\theta\rightarrow\theta+2\pi$ as a source
of the pion, as we did in our previous work Itou:2023img . We generate the
ground states for $2\pi\leq\theta\leq 3\pi$ and measure the one-point function
$\Braket{\pi(x)}$ of the pion (41). The result is shown in the left panel of
figure 19, which has similar behavior to the results with the twisted mass in
figure 6. However, we find $\Braket{\pi(x)}\approx 0$ at $\theta=3\pi$
unexpectedly. The reason would be that the SPT state with the boundary charge
gets higher energy than the trivial state due to the finite size effect around
the critical point $\theta=3\pi$. Although we cannot compare the result at
$\theta=3\pi$ with the analytic calculation as in section 6, we can still
obtain the pion mass for $2\pi\leq\theta<3\pi$. We compute the effective mass
(48) and fit the results by (47) in the range $10\leq x\leq 25$. The results
are plotted in the right panel of figure 19, and the fitting parameters are
summarized in table 7.

Figure 19:  (Left) The one-point function $\log|\Braket{\pi(x)}|$ of the pion is plotted against the distance $x$ from the boundary for $\varepsilon=10^{-10}$. (Right) The effective mass calculated from the one-point function in the left panel is plotted against $x$. The results of fitting by (47) are shown as well. The fitting range is between the vertical dashed lines. $\theta/2\pi$ | $M_{\pi}$ | $\Delta M$ | $C$
---|---|---|---
1.0 | 0.431167(3) | 0.4178(6) | -0.631(3)
1.1 | 0.418404(2) | 0.3981(5) | -0.772(3)
1.2 | 0.379502(5) | 0.3501(8) | -0.960(6)
1.3 | 0.31225(1) | 0.272(1) | -1.187(9)
1.4 | 0.2092(1) | 0.167(2) | -1.36(1)
Table 7:  The fitting results of the effective mass computed from the one-
point function $\Braket{\pi(x)}$ of the pion. The errors of these values come
from the fitting error.

Let us summarize the results of the one-point-function scheme with the naive
open boundary. The $\theta$-dependent masses of the pion and sigma meson are
shown in figure 20. The results are consistent with those obtained on the
lattice equipped with the wings in figure 7. So, they also agree with the
analytic results by the bosonization. In the one-point-function scheme, we
have many choices of the source of the mesons. In this work, we tried four
cases: the wings with the large fermion mass, the isospin-breaking chiral
rotation, the naive open boundary, and the edge mode of the SPT state. If we
can find an appropriate source for the target particle, the one-point-function
scheme should widely apply to other cases.

Figure 20:  The masses of the pion and sigma meson obtained by the one-point-
function scheme with the naive open boundary are plotted against
$\theta/2\pi$. The eta-meson mass at $\theta=0$ is also plotted for reference.
The gray dashed curve denotes the calculation by the bosonized model,
$M_{\pi}(0)|\cos(\theta/2)|^{2/3}$, where the overall coefficient $M_{\pi}(0)$
is determined by the numerical result at $\theta=0$. The gray dash-dot curve
denotes $\sqrt{3}M_{\pi}(0)|\cos(\theta/2)|^{2/3}$.

##  Appendix E Fate of the eta meson

At $\theta=0$, the eta meson is a stable particle protected by the $G$-parity.
However, the $G$-parity is explicitly broken for $\theta\neq 0$. Then the eta
meson is no longer stable but decays into the pions or mixes with the sigma
meson. In terms of the bosonized model, the eta meson has the interaction (12)
with the other mesons of the coupling $\sim\sin(\theta/2)$, as discussed in
section 2.2. Indeed, the eigenstates of the Hamiltonian corresponding to the
eta meson can be found at $\theta=0$, but such states disappear as we increase
$\theta$ in the dispersion-relation scheme. We observed that the eta-meson-
like state and the sigma meson are mixed in the spectrum depending on
$\theta$. In this appendix, we show the one-point and two-point correlation
functions of the eta meson and discuss how they behave when the eta meson is
unstable. We also show the values of the effective mass $M_{\eta}$ as a
reference when we naively apply the estimation method discussed in the main
text, and they should not be taken seriously as it is well-established only
for stable particles.

###  E.1 One-point function of the eta meson

First, we look into the one-point function. The boundary condition violating
the $G$-parity can be a source of the eta meson. Here we adopt the lattice
setup with the wings region assigned a large fermion mass
$m_{\mathrm{wings}}=m_{0}=10$ as in section 5.2. We measure the one-point
function $\Braket{\eta(x)}$ of the eta meson (41) for the ground states, using
the result in section 5.2.1 to determine the matrix $R_{+}$. The results of
$\Braket{\eta(x)}$ are shown in the left panel of figure 21. Unlike the pion
and sigma meson, the behavior of $\Braket{\eta(x)}$ is no longer a simple
exponential decay $\sim e^{-Mx}$ for $\theta\neq 0$, and it is not even
monotonically decreasing resulting in the cusp of the plot. The drastic change
of the behavior depending on $\theta$ should be caused by the decay or mixing.
Note that the eta-meson mass $M_{\eta}\sim 0.9$ is larger than twice the pion
mass $M_{\pi}\lesssim 0.43$, and thus the decay is allowed in the current
setup.

As long as the effect of decay is small, we can still measure the mass of the
unstable eta meson by examining the one-point function. Let us look into the
effective mass. The effective mass computed from $\Braket{\eta(x)}$ is shown
in the right panel of figure 21. As $\theta$ increases, the plateau region of
the effective mass shrinks. We fit the effective mass in the short range,
$3\leq x\leq 8$, by (47) for $\theta/2\pi=0.0,0.1,0.2$. The fitting results
are plotted in the right panel of figure 21, and the fitting parameters are
summarized in table 8. The $\theta$-dependence of the eta-meson mass at short
distances is relatively small compared with the masses of the pion and sigma
meson. At long distances, $x\gtrsim 25$, the effective mass takes lower values
of about 0.3 to 0.5 for $\theta/2\pi=0.3$ and $0.4$. This behavior should be
explained by the lighter states which the eta meson decays into. It is
interesting to check whether a similar behavior is observed in analytic
calculations assuming the effective model of the decaying eta meson.

Figure 21:  (Left) The one-point function $\log|\Braket{\eta(x)-\eta(L/2)}|$ of the eta meson is plotted against the distance $x$ from the boundary for $\varepsilon=10^{-10}$. The value at $x=L/2$ is subtracted from $\Braket{\eta(x)}$ to eliminate the constant shift in the bulk. The shaded region indicates the wings. (Right) The effective mass calculated from the one-point function in the left panel is plotted against $x$. The results of fitting by (47) are shown as well. The fitting range is between the vertical dashed lines. $\theta/2\pi$ | $M_{\eta}$ | $\Delta M$ | $C$
---|---|---|---
0.0 | 0.9022(1) | 0.787(7) | 0.029(1)
0.1 | 0.90319(5) | 0.814(3) | 0.0239(5)
0.2 | 0.89501(8) | 0.731(3) | 0.0348(8)
Table 8:  The fitting results of the effective mass at short distances
computed from the one-point function $\Braket{\eta(x)}$ of the eta meson. The
errors of these values come from the fitting error.

###  E.2 Correlation function of the eta meson

Next, let us focus on the spatial correlation function. The connected
correlation function of the eta meson,

$C_{\eta}(r)=\Braket{\eta(x)\eta(y)}_{c},$ (77)

is obtained as the eigenvalue of the correlation matrix (35) together with the
sigma meson. The results are shown in figure 22. In this plot, the behavior of
the correlation function for $\theta/2\pi\gtrsim 0.4$ is qualitatively
different from the other cases for the small $\theta$. The gradient of
$\log\Braket{\eta(x)\eta(y)}_{c}$ suddenly changes at long distances $r\gtrsim
10$ whereas it is less sensitive to $\theta$ at short distances $r\lesssim
10$.

Figure 22:  The correlation function of the eta meson
$\log|\Braket{\eta(x)\eta(y)}_{c}|$ is plotted against the distance $r=|x-y|$.
The cutoff parameter in DMRG is $\varepsilon=10^{-16}$.

Correspondingly, the behavior of the effective mass changes depending on the
distance $r$. The effective mass computed from the correlation function of the
eta meson is shown in figure 23 for $\theta/2\pi=0.0,0.2,0.4$. For the large
$\theta$ such as $\theta/2\pi=0.4$ in the right panel of figure 23, we find
that $M_{\eta,\mathrm{eff}}(r)$ behaves as $\alpha/r+M$ at both short and long
distances but with different $\alpha$ and $M$, which should result from the
decay and mixing. At short distances, the correlation function has information
on the unstable eta meson before decaying. At long distances, the correlation
function behaves like the lighter mesons, namely the pion after decaying or
the mixed sigma meson. From this point of view, we focus on the short-distance
region to obtain the mass of the unstable eta meson when $\theta$ is large. We
fit the effective mass $M_{\eta,\mathrm{eff}}(r)$ by $\alpha/r+M$ for
$\varepsilon=10^{-16}$. Here we evaluate the systematic errors from the
uncertainty of the fitting range by changing the range inside a maximum
fitting region, as in appendix C. The maximum region is $0.075\leq 1/r\leq
0.15$ for $\theta/2\pi=0.0,\cdots,0.3$. To see the short-distance behavior,
the maximum region is changed to $0.125\leq 1/r\leq 0.25$ for
$\theta/2\pi=0.4$ and $0.15\leq 1/r\leq 0.35$ for $\theta/2\pi=0.5$. The
fitting results are depicted by the purple shadows in figure 23, and the
results of the fitting parameter are summarized in table 9. The mass of the
unstable eta meson is less sensitive to $\theta$ than the other mesons.

Figure 23:  The effective masses $M_{\eta,\mathrm{eff}}(r)$ of the eta meson at $\theta/2\pi=0.0$ (left), $0.2$ (center), and $0.4$ (right) are plotted against $1/r$ for various values of $\varepsilon$. The data points for $\varepsilon=10^{-16}$ are fitted by $\alpha/r+M$. The fitting range is between the vertical dashed lines. The fitting result is depicted by the shaded band with the systematic error. For $\theta/2\pi=0.4$ we try two choices of the fitting range. The result at the short distance is used to obtain the mass of the unstable eta meson. $\theta/2\pi$ | $M_{\eta}$ | $\alpha$
---|---|---
0.0 | 0.899(3) | 0.51(3)
0.1 | 0.894(2) | 0.56(2)
0.2 | 0.872(4) | 0.68(4)
0.3 | 0.815(4) | 0.94(3)
0.4 | 0.720(3) | 1.30(2)
0.5 | 0.57(3) | 1.9(1)
Table 9:  The fitting results of the effective mass $M_{\eta,\mathrm{eff}}(r)$
computed from the correlation function of the eta meson. The errors of these
values come from the systematic error of fitting range.

Next, we investigate the effective mass at long distances, where the
contributions of the light mesons are expected to be significant. To estimate
the mass of the meson after decaying, we fit the effective mass
$M_{\eta,\mathrm{eff}}(r)$ by $\alpha/r+M$ for $\theta/2\pi=0.4$ and $0.5$
with different fitting ranges than before. For $\theta/2\pi=0.4$, we set the
maximum fitting region $0.05\leq 1/r\leq 0.08$ and obtain the results
$M=0.34(3)$ and $\alpha=5.7(5)$. The fitting result is depicted by the dark
yellow shadow in the right panel of figure 23. The obtained mass $M=0.34(3)$
is consistent with the sigma-meson mass at $\theta/2\pi=0.4$, for example,
$M_{\sigma}=0.334(5)$ by the correlation-function scheme in table 5. For
$\theta/2\pi=0.5$, we set the maximum region $0.06\leq 1/r\leq 0.125$ and
obtain $M=-0.02(2)$ and $\alpha=6.1(3)$. In this case, the mass is consistent
with zero, which agrees with the CFT-like behavior of the pion and sigma meson
at $\theta=\pi$. It is again interesting to compare these results with
analytic calculation assuming an effective model.

## References

  * (1) Flavour Lattice Averaging Group (FLAG) collaboration, _FLAG Review 2021_ , _Eur. Phys. J. C_ 82 (2022) 869 [2111.09849].
  * (2) K. Nagata, _Finite-density lattice QCD and sign problem: Current status and open problems_ , _Prog. Part. Nucl. Phys._ 127 (2022) 103991 [2108.12423].
  * (3) J.B. Kogut and L. Susskind, _Hamiltonian Formulation of Wilson’s Lattice Gauge Theories_ , _Phys. Rev. D_ 11 (1975) 395.
  * (4) E. Itou, A. Matsumoto and Y. Tanizaki, _Calculating composite-particle spectra in Hamiltonian formalism and demonstration in 2-flavor QED 1+1d_, _JHEP_ 11 (2023) 231 [2307.16655].
  * (5) S.R. White, _Density matrix formulation for quantum renormalization groups_ , _Phys. Rev. Lett._ 69 (1992) 2863.
  * (6) S.R. White, _Density-matrix algorithms for quantum renormalization groups_ , _Phys. Rev. B_ 48 (1993) 10345.
  * (7) U. Schollwöck, _The density-matrix renormalization group_ , _Reviews of Modern Physics_ 77 (2005) 259.
  * (8) U. Schollwöck, _The density-matrix renormalization group in the age of matrix product states_ , _Annals of Physics_ 326 (2011) 96.
  * (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) N.S. Manton, _The Schwinger Model and Its Axial Anomaly_ , _Annals Phys._ 159 (1985) 220.
  * (14) J.E. Hetrick and Y. Hosotani, _QED ON A CIRCLE_ , _Phys. Rev._ D38 (1988) 2621.
  * (15) C. Jayewardena, _SCHWINGER MODEL ON S(2)_ , _Helv. Phys. Acta_ 61 (1988) 636.
  * (16) I. Sachs and A. Wipf, _Finite temperature Schwinger model_ , _Helv. Phys. Acta_ 65 (1992) 652 [1005.1822].
  * (17) C. Adam, _Instantons and vacuum expectation values in the Schwinger model_ , _Z. Phys._ C63 (1994) 169.
  * (18) J.E. Hetrick, Y. Hosotani and S. Iso, _The Massive multi - flavor Schwinger model_ , _Phys. Lett._ B350 (1995) 92 [hep-th/9502113].
  * (19) R. Narayanan, _QED at a finite chemical potential_ , _Phys. Rev._ D86 (2012) 087701 [1206.1489].
  * (20) R. Narayanan, _Two flavor massless Schwinger model on a torus at a finite chemical potential_ , _Phys. Rev._ D86 (2012) 125008 [1210.3072].
  * (21) 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].
  * (22) Y. Tanizaki and M. Tachibana, _Multi-flavor massless QED 2 at finite densities via Lefschetz thimbles_, _JHEP_ 02 (2017) 081 [1612.06529].
  * (23) 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].
  * (24) 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].
  * (25) 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].
  * (26) 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].
  * (27) B. Buyens, J. Haegeman, H. Verschelde, F. Verstraete and K. Van Acoleyen, _Confinement and string breaking for QED 2 in the Hamiltonian picture_, _Phys. Rev. X_ 6 (2016) 041040 [1509.00246].
  * (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) B. Buyens, S. Montangero, J. Haegeman, F. Verstraete and K. Van Acoleyen, _Finite-representation approximation of lattice gauge theories at the continuum limit with tensor networks_ , _Phys. Rev. D_ 95 (2017) 094509 [1702.08838].
  * (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) D.E. Kharzeev and Y. Kikuchi, _Real-time chiral dynamics from a digital quantum simulation_ , _Phys. Rev. Res._ 2 (2020) 023342 [2001.00698].
  * (34) 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].
  * (35) S. Thompson and G. Siopsis, _Quantum computation of phase transition in the massive Schwinger model_ , _Quantum Sci. Technol._ 7 (2022) 035001 [2110.13046].
  * (36) 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].
  * (37) 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].
  * (38) L. Nagano, A. Bapat and C.W. Bauer, _Quench dynamics of the Schwinger model via variational quantum algorithms_ , _Phys. Rev. D_ 108 (2023) 034501 [2302.10933].
  * (39) L. Funcke, K. Jansen and S. Kühn, _Exploring the CP-violating Dashen phase in the Schwinger model with tensor networks_ , _Phys. Rev. D_ 108 (2023) 014504 [2303.03799].
  * (40) T. Angelides, L. Funcke, K. Jansen and S. Kühn, _Computing the mass shift of Wilson and staggered fermions in the lattice Schwinger model with matrix product states_ , _Phys. Rev. D_ 108 (2023) 014516 [2303.11016].
  * (41) K. Ikeda, D.E. Kharzeev, R. Meyer and S. Shi, _Detecting the critical point through entanglement in the Schwinger model_ , _Phys. Rev. D_ 108 (2023) L091501 [2305.00996].
  * (42) 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_ , _Phys. Rev. Lett._ 132 (2024) 031603 [2305.04437].
  * (43) J.W. Pedersen, E. Itou, R.-Y. Sun and S. Yunoki, _Quantum Simulation of Finite Temperature Schwinger Model via Quantum Imaginary Time Evolution_ , _PoS_ LATTICE2023 (2024) 220 [2311.11616].
  * (44) P. Schmoll, J. Naumann, A. Nietner, J. Eisert and S. Sotiriadis, _Hamiltonian truncation tensor networks for quantum field theories_ , 2312.12506.
  * (45) T. Angelides, P. Naredi, A. Crippa, K. Jansen, S. Kühn, I. Tavernelli et al., _First-Order Phase Transition of the Schwinger Model with a Quantum Computer_ , 2312.12831.
  * (46) D. Ghim and M. Honda, _Digital Quantum Simulation for Spectroscopy of Schwinger Model_ , _PoS_ LATTICE2023 (2024) 213 [2404.14788].
  * (47) P.P. Popov, V. Kasper, M. Lewenstein, E. Zohar, P. Stornati and P. Hauke, _Non-perturbative signatures of fractons in the twisted multi-flavor Schwinger Model_ , 2405.00745.
  * (48) O. Kaikov, T. Saporiti, V. Sazonov and M. Tamaazousti, _Phase Diagram of the Schwinger Model by Adiabatic Preparation of States on a Quantum Simulator_ , 2407.09224.
  * (49) 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].
  * (50) Y. Shimizu and Y. Kuramashi, _Grassmann tensor renormalization group approach to one-flavor lattice Schwinger model_ , _Phys. Rev. D_ 90 (2014) 014508 [1403.0642].
  * (51) Y. Shimizu and Y. Kuramashi, _Critical behavior of the lattice Schwinger model with a topological term at $\theta=\pi$ using the Grassmann tensor renormalization group_, _Phys. Rev. D_ 90 (2014) 074503 [1408.0897].
  * (52) N. Butt, S. Catterall, Y. Meurice, R. Sakai and J. Unmuth-Yockey, _Tensor network formulation of the massless Schwinger model with staggered fermions_ , _Phys. Rev. D_ 101 (2020) 094509 [1911.01285].
  * (53) A. Yosprakob, J. Nishimura and K. Okunishi, _A new technique to incorporate multiple fermion flavors in tensor renormalization group method for lattice gauge theories_ , _JHEP_ 11 (2023) 187 [2309.01422].
  * (54) F.I. Az-zahra, S. Takeda and T. Yamazaki, _Spectroscopy by Tensor Renormalization Group Method_ , 2404.15666.
  * (55) 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].
  * (56) 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].
  * (57) 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].
  * (58) T. Sulejmanpasic and C. Gattringer, _Abelian gauge theories on the lattice: $\theta$-terms and compact gauge theory with(out) monopoles_, _Nucl. Phys._ B943 (2019) 114616 [1901.02637].
  * (59) T. Sulejmanpasic, D. Göschl and C. Gattringer, _First-Principles Simulations of 1+1D Quantum Field Theories at $\theta=\pi$ and Spin Chains_, _Phys. Rev. Lett._ 125 (2020) 201602 [2007.06323].
  * (60) H. Ohata, _Monte Carlo study of Schwinger model without the sign problem_ , _JHEP_ 12 (2023) 007 [2303.05481].
  * (61) H. Ohata, _Phase diagram near the quantum critical point in Schwinger model at $\theta$ = $\pi$: analogy with quantum Ising chain_, _PTEP_ 2024 (2024) 013B02 [2311.04738].
  * (62) S.R. Coleman, _More About the Massive Schwinger Model_ , _Annals Phys._ 101 (1976) 239.
  * (63) S.R. Coleman, _The Quantum Sine-Gordon Equation as the Massive Thirring Model_ , _Phys. Rev. D_ 11 (1975) 2088.
  * (64) 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.
  * (65) I. Affleck and F.D.M. Haldane, _Critical Theory of Quantum Spin Chains_ , _Phys. Rev._ B36 (1987) 5291.
  * (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) M.L. Wall and L.D. Carr, _Out-of-equilibrium dynamics with matrix product states_ , _New Journal of Physics_ 14 (2012) 125015.
  * (69) M. Fishman, S.R. White and E.M. Stoudenmire, _The ITensor Software Library for Tensor Network Calculations_ , _SciPost Phys. Codebases_ (2022) 4.