Calculating composite-particle spectra in Hamiltonian formalism and demonstration in 2-flavor QED1+1d

We consider three distinct methods to compute the mass spectrum of gauge theories in the Hamiltonian formalism: (1) correlation-function scheme, (2) one-point-function scheme, and (3) dispersion-relation scheme. The first one examines spatial correlation functions as we do in the conventional Euclidean Monte Carlo simulations. The second one uses the boundary effect to efficiently compute the mass spectrum. The third one constructs the excited states and fits their energy using the dispersion relation with selecting quantum numbers. Each method has its pros and cons, and we clarify such properties in their applications to the mass spectrum for the 2-flavor massive Schwinger model at m/g = 0.1 and θ = 0 using the density-matrix renormalization group (DMRG). We note that the multi-flavor Schwinger model at small mass m is a strongly coupled field theory even after the bosonizations, and thus it deserves to perform the first-principles numerical calculations. All these methods mostly agree and identify the stable particles, pions πa (JPG = 1−+), sigma meson σ (JPG = 0++), and eta meson η (JPG = 0−−). In particular, we find that the mass of σ meson is lighter than twice the pion mass, and thus σ is stable against the decay process, σ → ππ. This is consistent with the analytic prediction using the WKB approximation, and, remarkably, our numerical results are so close to the WKB-based formula between the pion and sigma-meson masses, Mσ/Mπ = 3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{3} $$\end{document}.

In recent years, numerical simulations of quantum field theories (QFTs) in the Hamiltonian formalism have attracted a lot of attention motivated by the rapid progress of quantum computing technology and also the developments of tensor network techniques.These methods rely on different disciplines from that of Monte Carlo simulations for the conventional lattice gauge theories, and thus they are expected to give complementary frameworks.One of the remarkable features is that these methods do not rely on importance sampling, and thus we may be able to circumvent the issue of sign problems.
With this motivation in mind, we investigate the methods to calculate the mass spectrum of gauge theories in Hamiltonian formalism.When studying strongly coupled QFTs, we often encounter situations where the fundamental degrees of freedom defining the theory do not appear in the low-energy spectrum.Quantum chromodynamics (QCD) is a notable, successful example of such phenomena: Quarks and gluons are confined inside the color-singlet hadrons, and it explains the physics of strong interaction.The Monte Carlo simulations nicely predict the hadron spectrum [1] and also the physics at finite temperature [2,3] in the sign-problem-free regions.Of course, we are currently very far away to reproduce such tremendous achievements of Monte Carlo simulations, and thus it is important to develop the counterparts of those calculational techniques in Hamiltonian formalism.
In this work, we consider three independent methods to compute the mass spectrum of the 2-flavor massive Schwinger model using the density-matrix renormalization group (DMRG): • correlation-function scheme • one-point-function scheme

• dispersion-relation scheme
The first one examines spatial correlation functions as we do in the conventional Euclidean Monte Carlo simulations.The second one uses the boundary effect for efficiently computing the mass spectrum, which is partly motivated by the applications of the Friedel oscillations [4,5].The third one constructs the excited states and fits their energy using the dispersion relation with selecting quantum numbers (see, e.g., Refs.[6][7][8] for similar analysis in spin systems).Each method has its pros and cons, and especially the last one is specific to the Hamiltonian formalism.Our purpose is to clarify their properties in the concrete studies of the 2-flavor massive Schwinger model.
The low-energy mass spectrum of the 2-flavor Schwinger model with a theta term has been studied analytically by using the bosonization technique [13,20].We note, however, that the low-energy effective theory is strongly coupled if 0 < m/g ≪ 1 even after bosonization, and the details of the prediction rely on some approximations that are not fully justifiable.It is physically nontrivial if those analytic predictions are reproduced by the first-principles numerical computations when we go into the details beyond qualitative aspects.
We performed the DMRG computations using the C++ library of ITensor [48] to obtain the mass spectrum at θ = 0 with the above three methods.We find that all three methods mostly agree and identify the stable particles, pions π a (J P G = 1 −+ ), sigma meson σ (J P G = 0 ++ ), and eta meson η (J P G = 0 −− ), where J denotes the isospin quantum number, P and G is the parity and G-parity, respectively.In particular, we observe that the mass of σ meson is lighter than twice the pion mass, and thus σ is stable against the decay process, σ → ππ.This implies that σ is a stable particle, not a ππ resonance, and this is a notable difference compared with 4d QCD.This is consistent with the analytic prediction based on the WKB approximation of the Abelian bosonized description.The WKB-based formula predicts the π and σ mass ratio is given by M σ /M π = √ 3, and our three distinct computations give roughly consistent results: 1 correlation-function one-point-function dispersion-relation M σ /M π 1. 68(2) 1.821 (6) 1.75(1) (1.1) Let us emphasize that this poses an interesting theoretical question on why the semiclassical prediction works so well even outside its valid regime.In this paper, we concentrate on the simulation at θ = 0 to confirm the validity of our method.In fact, it should be straightforward to extend our work to the θ ̸ = 0 region.This paper is organized as follows.In Section 2, we review the continuum 2-flavor Schwinger model and the bosonization analysis focusing on the mass spectrum.In Section 3, we introduce the lattice formulation of the Hamiltonian and define some observables.In Section 4, we briefly explain our simulation method, the DMRG algorithm and show our setup of the simulation.In Section 5, we present our simulation results for the three methods.Section 6 is devoted to the conclusion and discussion.Appendix A shows the explicit form of the Hamiltonian and the observables in the spin representation for DMRG.In Appendix B, we test the validity of the charge conjugation operator on the lattice in the 1-flavor Schwinger model.In Appendix C, we discuss the assignment of the flavor index for constructing MPS.In Appendix D, we investigate how the truncation of the bond dimension affects the correlation function in the massless 1-flavor Schwinger model.

Review of the 2-flavor Schwinger model
In this work, we study the 2-flavor Schwinger model, which is (1 + 1)-dimensional quantum electrodynamics (QED 1+1d ) with N f = 2 species of Dirac fermion.The Lagrangian density with the Minkowski metric η µν = diag(1, −1) is given by where is the field strength, g is the gauge coupling, and θ is the vacuum angle describing the background electric flux.The index f labels the flavor.We set the masses of the two fermions equal to m.

Global symmetry and composite operators
In the chiral limit (m = 0), the 2-flavor Schwinger model has the chiral symmetry and the G-parity symmetry, SU(2 and the chiral symmetry has an 't Hooft anomaly.The G-parity operation is the combination of the charge-conjugation with the π rotation of the SU(2) V , which will be discussed later.We note that the continuous chiral symmetry cannot be spontaneously broken due to the Coleman-Mermin-Wagner theorem, and the anomaly matching condition is satisfied by the SU(2) level-1 Wess-Zumino-Witten (SU(2) 1 WZW) conformal field theory.The SU(2) 1 WZW model is equivalent to the self-dual compact boson, and one can explicitly derive it from the massless 2-flavor Schwinger model with the Abelian bosonization [13].
The massive 2-flavor Schwinger model (2.1) no longer has the chiral symmetry, but it maintains the following symmetry, and we call SU(2) V /Z 2 as the isospin symmetry.(Z 2 ) G+L is the diagonal subgroup between the center of SU(2) L and the G-parity (Z 2 ) G .The Z 2 quotient of SU(2)/Z 2 means that the local operators always have the integer isospin quantum numbers since the gauge invariance requires that the local operator must consist of the same number of ψ and ψ.We define the isospin operators J a as conserved charges under this symmetry by where τ a represents Pauli matrices of the isospin space with a ∈ {x, y, z}.
When θ takes some special values, e.g.θ = 0, the theory enjoys the charge conjugation, with a suitable element of the Clifford algebra C. We note that this operation flips the sign of the θ angle, and thus this symmetry does not exist for generic values of θ.For general numbers of flavors, C acts on SU(N f )/Z N f as an outer automorphism, i.e. the symmetry group becomes does not have nontrivial outer automorphisms, and indeed C just gives the π rotation in the isospin space.Thus, it is convenient to introduce the G-parity [49], so that it commutes with the isospin operation and gives a well-defined eigenvalue ±1.
Moreover, the G-parity acts trivially on the SU(2) 1 WZW theory.Thus, if we find a particle with G = −1, we can immediately tell it remains massive in the chiral limit.
In this paper, we mainly focus on the following composite operators to discuss the meson spectrum: ) (2.9) We call them pion, sigma, and eta operators, respectively, obviously motivated by the meson spectrum of 4d QCD.We will often denote π 3 = π for simplicity.Here, we have specified their quantum numbers J P G , where J denotes the isospin, and P and G denote the parity and the G-parity at θ = 0, respectively.The (1 + 1)d QED is strongly coupled, and it turns out that the light particles correspond to these operators, and this feature is reminiscent of 4d QCD.Here, we would like to note that η has G = −1, and thus it remains massive in the chiral limit, which is analogous to the U (1) A problem [12,13,50].
Other mesons, π and σ have G = +1 and actually become massless in the chiral limit.The massless σ particle is an outcome of the absence of chiral symmetry breaking, unlike the 4d QCD case.

Phase structure
With m ̸ = 0, the system is gapped and has the unique ground state at generic values of θ.
In (1 + 1)d, there is no stable topologically-ordered state, and the unique gapped ground states are then classified as the symmetry-protected topological (SPT) states [51][52][53].This perspective provides a very powerful tool to understand the phase structure of the 2-flavor Schwinger model.Let us recall that the massive Schwinger model always has the isospin symmetry, SU(2) V /Z 2 .We can then calculate the partition function under the presence of the background gauge fields for the isospin symmetry.Compared with the SU(2) gauge field, the SU(2)/Z 2 gauge field has milder cocycle conditions, which is controlled by the Z 2 2-form gauge field w 2 in addition to the familiar 1-form gauge field.As a result, at generic values of θ, the partition function with the background gauge field is described by the low-energy effective topological action, with some k ∼ k + 2. We note that k is a discrete label that distinguishes the SPT states protected by SU(2) V /Z 2 , and thus it cannot be changed under the continuous change of coupling constants unless quantum phase transitions happen.We can prove that two partition functions at θ and θ + 2π are different in the presence of the background gauge fields.The anomalous relation can be summarized as [54] Z θ+2π = exp iπ w 2 Z θ . (2.11) The label k is changed as k → k + 1 as we change the θ angle from θ to θ + 2π, and there must be a phase transitions separating the k = 0, 1 ground states.It is somewhat customary to assign the k = 0 state for −π < θ < π and the k = 1 state for π < θ < 3π, while, precisely, this assignment depends on the UV-regularization scheme.We note that the whole story here is quite parallel to that of the anti-ferromagnetic Heisenberg chain or the (1 + 1)d CP 1 sigma model [55][56][57][58][59][60][61][62] (except for the properties at θ = π [13,38]).The distinction between the states at θ = 0 and θ = 2π becomes more vivid when we take the open boundary condition.Turning on non-zero θ corresponds to introducing a background electric field with a constant magnitude θ/2π.When we increase θ beyond π, the background field becomes larger than 1/2.Then the Dirac fermions with charges ±1 are excited at the boundaries to cancel the background field as much as possible.As a consequence, these boundary states have isospin 1/2, which is the projective representation of SU(2) V /Z 2 .This is nothing but the signature of the nontrivial SPT state.If the system size is large enough, the interaction between the boundary states is exponentially suppressed, so the two independent degrees of freedom with isospin 1/2 yield 2 × 2 degeneracy of the ground state.We will use this boundary excitation as a source of iso-triplet particles to determine the mass of the pion from the one-point function in Section 5.2.Note that if we increase θ further beyond 3π, the boundary excitations become bound states of Dirac fermions with isospin 1, which can be completely screened by the gauge-invariant particles inside the bulk, and the ground state would be unique again with the open boundary condition.

Mass spectrum
In this subsection, we are going to give a relatively detailed review on analytic predictions about the mass spectrum.There is already a huge effort for the analytic studies of the multiflavor Schwinger model, and thus the reader may wonder how one could obtain something new with numerical studies.We would like to clarify what kinds of approximations were used in the previous studies and give the justification of this work on physical aspects.
There are two exactly-solvable limits of the multi-flavor Schwinger model: • In the heavy fermion limit m → ∞, the model becomes the pure U (1) gauge theory.
• In the chiral limit m = 0, the model becomes the SU (N f ) 1 WZW conformal field theory and one massive free boson, and they are completely decoupled.
It is then natural to consider perturbations from these limits in order to investigate the cases of general mass m.When m ≫ g, one can perform systematic perturbations to study the spectrum.In the opposite case 0 < m ≪ g, however, the systematic perturbation works only for the 1-flavor case, and further approximations are necessary for N f ≥ 2. By applying the Abelian bosonization for N f = 2, the fermions are mapped to the 2π-periodic scalar fields ϕ 1 , ϕ 2 .The Lagrangian (2.1) is then completely equivalent to where C = e γ /(2π) is a numerical constant, and N ρ denotes the normal ordering for the contraction with a free field propagator of mass ρ [63]. 2 We can rigorously integrate out the gauge fields as the Lagrangian is quadratic in terms of F 01 .Changing the basis of bosons as ϕ 1,2 = √ 2πη − θ 2 ± φ, the effective Lagrangian becomes where µ 2 = 2g 2 /π.We have the Z 2 symmetry acting only on η when θ = 0, and this is the G-parity.When m = 0, the massive η and the massless φ decouple, as advocated above.
Let us now turn on the small mass, 0 < m ≪ g.The η particle has the mass µ + O(m), but it is hard to compute the O(m) correction due to the potential infrared divergence in the loop diagrams with the φ fields [13].Instead, we integrate out η at the tree level to discuss the physics of π and σ mesons, which gives N ρ cos √ 2πη − θ 2 = µ ρ cos θ 2 for the free massive η with mass µ.The effective theory for φ becomes the sine-Gordon model, The isospin SU(2) V /Z 2 symmetry is not manifest at all in the Lagrangian, but it secretly exists quantum mechanically.In particular, the z-component isospin current is given by J µ z = 1 2π ε µν ∂ ν φ, and its charge J z = dx 1 2π ∂ x φ counts the winding number of φ ∼ φ + 2π.We would like to emphasize that the effective theory (2.14) is strongly coupled, and we cannot solve it in the ordinary perturbation for small but nonzero m.What is actually done in the previous literature is the optimized perturbation; we optimize the renormalization scale ρ so that the coefficient of the cos(φ) potential becomes O(ρ 2 ), and we get This is identified as the mass gap caused by the mass perturbation, and this formula gives the θ-dependence of the lightest meson mass, i.e.M π .
The spectrum of the sine-Gordon model was studied by using WKB approximation [64].Introducing an extra parameter controlling the kinetic term as 1 4πβ 2 (∂φ) 2 , the quantum scaling dimension of cos φ becomes ∆ = β 2 /2, so the semiclassical approximation is valid when β 2 → 0. The model has the soliton and antisoliton, and let us denote their mass as M SG .Then, Dashen et al. [64] predicted the masses of soliton-antisoliton bound states as ) . Even though it is subtle if the WKB works at the self-dual point β 2 = 1, Coleman got an intriguing observation using this semiclassical formula [13].
The nontrivial check for its validity is the recovery of the SU(2) V /Z 2 symmetry at β 2 = 1.Substituting n = 1 with β 2 = 1 into (2.16),we get (2.17) This shows that the lightest soliton-antisoliton bound state has the same mass as the soliton or antisoliton itself.The soliton and antisoliton have J z = ±1, and the soliton-antisoliton bound state has J z = 0.The G-parity does not act on φ, so all these states have G = +1.Thus, these three states form the isospin triplet J P G = 1 −+ of mass M SG , which is identified as the pion in the Schwinger model, and then . The mass of the second soliton-antisoliton bound state is given by This state has J z = 0, and there is no other state with the same mass.We then identify it as the σ meson in the Schwinger model with J P G = 0 ++ .Thus, the semiclassical method predicts that the masses of pion and sigma meson satisfy Importantly, M σ < 2M π .Unlike the 4d QCD, σ is a stable particle, not a resonance, because the decay σ → ππ is energetically prohibited.As we discussed above, the low-energy mass spectra can be estimated by bosonization.However, it relies on the optimized perturbation and also on the semiclassical method, and these analyses are not necessarily fully justified.It is still difficult to compute the exact m-dependence or to find other states with higher energies than µ.Thus, it is worth studying the mass spectrum of the 2-flavor model by first-principles numerical methods.

Lattice formulation of the 2-flavor Schwinger model
In this section, we explain the Hamiltonian formalism of the 2-flavor Schwinger model and its lattice regularization as a generalization to N f = 2 from previous research [32][33][34][35].We also define various local and global observables used in the analysis.

Hamiltonian
First, we introduce the continuum Hamiltonian of the N f -flavor Schwinger model.By introducing a conjugate momentum Π = 1 g 2 ∂ 0 A 1 + θ 2π , the Hamiltonian is given by In Hamiltonian formalism, the physical Hilbert space is constrained by the Gauss law condition, The electric field corresponds to E := Ȧ1 = g 2 (Π−θ/2π).Thus, the theta angle θ plays the role of the background electric field.In the periodic boundary condition, the Hamiltonians at θ and θ + 2π are unitary equivalent, realizes the 2π periodicity of θ.
Next, we consider the lattice regularization of the Hamiltonian.Here we employ the staggered fermion to define fermions on the lattice [65,66].The staggered fermions χ f,n with the lattice spacing a represents the discretization of the two-component Dirac fermions ψ f (x) with the lattice spacing 2a.The single-component fermions χ f,n at the site 3) The number of staggered fermions for each flavor is equal to N , thus at each site, there are two staggered fermions.In this work, we set N to be an even number.
The gauge field is encoded to U(1) variables U n ∼ exp −iaA 1 (x) , defined on the link between the n-th and (n + 1)-th sites, and the conjugate momentum is replaced by L n ∼ −Π(x), defined on the n-th site.The canonical commutation relations are given by ) Note that the roles of the staggered fermion operators depend on the site n: χ † f,n : creation op. of particle n : even annihilation op. of anti-particle n : odd χ f,n : annihilation op. of particle n : even creation op. of anti-particle n : odd Thus, the operator χ † f,n χ f,n counts the number of particles on the even sites, whereas χ f,n χ † f,n counts the number of anti-particles on the odd sites.Considering that the particle has an electric charge of +1 and the anti-particle has −1, the charge density operator at the site n is given by In this work, we choose the open boundary condition in order to eliminate the bosonic degrees of freedom having an infinite dimensional Hilbert space.The Gauss law (3.2) is also discretized as where the left-hand side corresponds to the divergence of the electric field and the righthand side is the charge density (3.9).We set the explicit form of the (1+1)d gamma matrices γ 0 = σ 3 , γ 1 = iσ 2 and γ 5 = γ 0 γ 1 = σ 1 .Using the operators χ f,n , U n , and L n introduced above, the lattice Hamiltonian is given by [37,38] where J = g 2 a/2 and w = 1/2a.Here we replace the mass m of the continuum theory by in the lattice Hamiltonian, following the recent proposal [67] for eliminating O(a) correction.In the continuum theory, the chiral limit m = 0 has the continuous chiral symmetry, With the above replacement, the lattice theory at m = 0 maintains the discrete chiral symmetry Z 2 ⊂ (Z N f ) L for even N f , and this is the point protected by the remnant of the chiral symmetry.By adding up the lattice Gauss law equation (3.10) from the boundary n = 0 to the site n, we find that L n can be replaced by where we set L −1 = 0. Furthermore, we can set U n = 1 since the degrees of freedom of U n can be absorbed by the U(1) phase of χ n .Then the lattice Hamiltonian is written only by the fermions as H = H J + H w + H m , where the gauge part H J is given by and the kinetic term H w and the mass term H m of the fermions are

Map to the spin system
Now, we map the Hamiltonian written by the staggered fermions to the spin Hamiltonian.Such a spin Hamiltonian formalism is useful to apply tensor network methods and quantum computations.
The N f × N degrees of freedom of the staggered fermion χ f,n can be described by the same number of spin-1/2 degrees of freedom.The Hilbert space of such a spin system is given by where H f,n is the local Hilbert space of the single spin-1/2 state.A general state |Ψ⟩ in this Hilbert space can be described by a superposition of all possible spin configurations s, The spin up |↑⟩ and down |↓⟩ state are the eigenstates of the Pauli matrix σ z with the eigenvalues +1 and −1, respectively.The map to the spin system can be achieved by the so-called Jordan-Wigner transformation.The fermion operators χ f,n for the two flavors f = 1, 2 are represented by spin operators as follows: where we define The Pauli matrices σ a f,n (a = x, y, z) act on the spin |s f,n ⟩ f,n at the site n of the flavor f .They do not commute only if they are on the same site of the same flavor, so that We can check that the canonical anti-commutation relations (3.4) and (3.5) are satisfied thanks to the properties of the Pauli matrices.Note that this is not a unique way of translation to the spin system which realizes the anti-commutation relations.Different transformations give different representations of the original Hamiltonian.We choose this transformation since various local operators can be constructed by only a few numbers of the Pauli matrices.The spin representation of the Hamiltonian and the observables defined above are summarized in Appendix A.

Local observables
Let us consider the meson operators (2.7) -(2.9) on the lattice.Based on the continuum descriptions, it is natural to define the lattice version of these operators as follows: Here S f,n and P S f,n are the scalar and pseudo-scalar operators for the flavor f = 1, 2 on the lattice, respectively.In order to obtain their explicit form, we rewrite the scalar condensate ( ψψ) f by the staggered fermion, so that Similarly, the pseudo-scalar condensate ( ψγ 5 ψ) f is given by These operators have a site-by-site fluctuation due to the staggered fermion.Here we define the lattice scalar condensate operator by the two-site average of (3.27), namely The lattice pseudo-scalar condensate operator is also defined by the two-site average of (3.28) with a factor −i, Note that both of S f (n) and P S f (n) are composed of the staggered fermions at the three sites n and n ± 1.

Global observables
We will define the quantum number (J z , J 2 , C, and P ) and momentum operators, which will be useful to distinguish the eigenstates of the Hamiltonian.These operators can be described by some global observables, which act on the whole lattice.
First of all, let us focus on the isospin operator (2.4).We define the lattice version in terms of the staggered fermion.The isospin J z operator counts the number of particles of each flavor with the factor ±1/2 on even sites and the number of anti-particles with the opposite sign on odd sites.Thus, it can be realized by It is convenient to define the isospin J ± operators by Based on the role of fermion operators (3.7) and (3.8), J + operator is given by which transforms f = 2 particle to f = 1 particle on even sites and f = 1 anti-particle to f = 2 anti-particle on odd sites.Similarly, J − operator is given by which transforms f = 1 particle to f = 2 particle on even sites and f = 2 anti-particle to f = 1 anti-particle on odd sites.Then the Casimir operator J 2 can also be defined as the combination of the operators above by Second, we will consider the charge conjugation and parity operators.For this purpose, let us discuss the description of the particle and anti-particle as a spin state.Applying the Jordan-Winger transformation, the spin representation of the charge density operator (3.9) is given by This operator counts the number of particles with +1 on even sites and the number of anti-particles with −1 on odd sites.We can confirm that the particle is described by the spin-up state |↑⟩ on even sites by taking the expectation value Similarly, we find that the anti-particle is described by the spin-down state |↓⟩ on odd sites as Based on this fact, charge conjugation, namely the exchange of particles and antiparticles can be performed by the exchange of even sites and odd sites.In addition, the spin-up state should be replaced by the spin-down state, and vice versa.These operations can be realized by the 1-site translation of the lattice and the multiplication of σ x operators.Thus, we define the charge conjugation operator by [25] where the swap operator is given by using the Pauli matrices.As the name suggests, the swap operator exchanges the state |s⟩ f,j and |s ′ ⟩ f,k , namely The product of the swap operators in (3.39) realizes the 1-site translation.The charge conjugation defined in this way satisfies When we take the periodic boundary condition, C 2 = 1 is achieved in the continuum limit, but this is not the case for the open boundary condition.Moreover, the Hamiltonian does not commute with C due to the presence of the boundaries, and we will actually see the eigenstates of the Hamiltonian give | ⟨C⟩ | < 1.Therefore, C does not give a good quantum number when we take the staggered-fermion regularization with the open boundary condition.In this study, following the observation of Ref. [25], we assume that the sign of Re ⟨C⟩ remembers the original sign of C for each eigenstate.We discuss this prescription in detail in Appendix B. Next, we define the parity operator.The parity transformation x → −x can be achieved by flipping the order of the lattice sites.The site n ∈ {0, 1, • • • , N − 1} is mapped to the site n ′ = N − 1 − n.However, this operation also exchanges particles and antiparticles since the roles of even sites and odd sites are exchanged when N is even.Thus, an additional operation of 1-site translation is necessary to fix it.We define the parity operator by where the products of the swap operators perform the reversal n → n ′ and the 1-site translation. 4The additional factors of σ z come from the shift of the staggered phase, which corresponds to γ 0 in the parity transformation of the Dirac fermion ψ(x) → γ 0 ψ(−x).As we mentioned for the C operator, the P operator in the open boundary condition does not commute with the Hamiltonian as it contains the 1-unit lattice translation.Therefore, we take the same prescription to determine the parity quantum number for each state as in the case of C. Finally, the other important quantity is a total momentum operator, which can be used to identify the momentum excitation [25].We start with the continuum description of the gauge invariant operator, which commutes with the continuum Hamiltonian (3.1) under the periodic boundary condition using the Gauss-law constraint (3.2).In our case with the open boundary condition, it does not commute with the Hamiltonian since the translational symmetry is explicitly broken.Thus, the expectation value ⟨K⟩ is no longer the quantum number in the strictest sense.However, we will see that the operator is still useful as an approximate one to investigate the mass spectrum of the model.Let us consider its lattice version.Here we set A 1 (x) = 0 since we fix the gauge U n = 1 in our setup.The combination ψ † f ∂ x ψ f of the Dirac fermion corresponds to in terms of the staggered fermion.There is another possible combination given by the integral by parts ignoring boundary term.Then we define the total momentum on the lattice as a Hermitian operator by taking symmetric combination

Calculation method and the simulation setup
We employ the density-matrix renormalization group (DMRG) [68][69][70][71] to study the spin Hamiltonian of the 2-flavor Schwinger model after the Jordan-Winger transformation, whose explicit form is given by (A.7).The DMRG is known as an efficient method to study (1+1)d gapped spin systems and has been developed mainly in the field of condensed matter physics.We utilized the C++ library of ITensor [48] to perform the tensor network calculation of this work.Let us briefly explain the basic idea of DMRG to obtain the ground state and excited states, and then we explain the details of parameter settings.

Quick review of DMRG
In the spin systems, any wave function |Ψ⟩ can be expressed as the form of the matrix product states (MPS), by repeating the singular-value decomposition (SVD).Here, i n = 1, 2 denotes the spin degrees of freedom at the n-th site, A (in) n denotes a D × D matrix, and this size D is called the bond dimension.The upper bound for the entanglement entropy of |Ψ⟩ is given by ln D. Therefore, the MPS gives a useful tool to study the many-body states with low entanglement entropies, such as the ground state of the (1 + 1)d gapped systems [72,73].For the 2-flavor Schwinger model, we arrange the site index n and the flavor index f on the 1d lattice with the single index to apply DMRG.The ordering of the indices is chosen so that the behavior of entanglement entropy is reproduced appropriately with a reasonable bond dimension.This point is discussed in Appendix C.
The DMRG is a variational algorithm based on the MPS.In each step of the algorithm, the matrices are updated to decrease the energy E = ⟨Ψ | H | Ψ⟩ as a cost function.In addition, we perform the low-rank approximation and thus the smaller singular values are discarded, which amount to an error ∆.We determine the bound dimension by setting the maximal bond dimension and also the cutoff parameter ε on the error so that ∆ ≤ ε.Smaller ε gives a better approximation, but it also requires a larger bond dimension and increases the computational costs.We can also effectively calculate the expectation values or correlation functions of local operators by rewriting those operators in the form of matrix product operators (MPOs) and then taking contractions with the ground state |Ψ⟩. 5e can use DMRG to obtain the low-energy excited states in a recursive way.Assume that we already find the energy eigenstates |Ψ ℓ ′ ⟩ with ℓ ′ = 0, 1, . . ., ℓ − 1 from below.Then, we apply the same technique to find the ℓ-th state |Ψ ℓ ⟩ by changing the Hamiltonian for the cost function as where W > 0 is a weight to impose the orthogonality.We can generate the excited states from the ground state to any level step by step.

Simulation setup
Let us explain our parameter setup when using the ITensor [48].The gauge coupling g has mass dimension 1 in 1 + 1d QED, and thus we can measure the energy scale in the unit of g by setting g = 1.In this work, we always set g = 1 and the fermion mass m = 0.1, so the photon mass is µ = 2 π ≃ 0.8.The lattice fermion mass (3.12) becomes m lat = 0.1 − a 4 .The theta angle is normally set to θ = 0, except when measuring the one-point function of the pion at θ = 2π.
For the correlation-function scheme and the one-point-function scheme, we use the lattice size of N = 160.The lattice spacing is set to a ≈ 0.25 so that the physical size is L = a(N − 1) = 39.8.The number of DMRG steps called the sweeps, is set to N sweep = 20.We generate the ground state for four different values of the cutoff parameter: ε = 10 −10 , 10 −12 , 10 −14 , and 10 −16 .To characterize the bond dimension of the MPS, we focus on the largest number of nonzero singular values, which we call the effective bond dimension, denoted as D eff .In our computations, we set the maximal bond dimension large enough so that D eff is solely controlled by the cutoff ε for the above physical setup.We observe D eff to be about 400, 800, 1600, and 2800 for the respective values of ε above.
For the dispersion-relation scheme, we generate many excited states up to ℓ = 23, which require a lot of computational cost.Therefore, we choose a smaller lattice size of L = 19.8 with N = 100 and a = 0.2.The excited states are generated with a cutoff of ε = 10 −10 and a weight parameter of W = 10.To achieve better convergence of higher states, we increase the number of sweeps to N sweep = 50.The bond dimension is about 500 for the ground state while it is at most 2300 for the excited states.
As an initial state of the DMRG, we choose the Néel state, which is a direct product of the spin-down states on even sites and the spin-up state on odd sites, Based on (3.37) and (3.38), the Néel states is regarded as a zero-particle state.We also impose the charge conservation condition during the DMRG, so that the MPS satisfies the condition Q = 0, where We note that the Gauss law with the usual open boundary on both sides requires Q = 0 on the physical states.

Simulation results
In this section, we explain our numerical results for the meson spectrum of the 2-flavor Schwinger model at θ = 0. We apply three distinct methods in our computations of the meson spectrum: • the correlation-function scheme • the one-point-function scheme • the dispersion-relation scheme Each method has its own pros and cons, and we are going to discuss them.We will see that all these schemes give consistent results.

Correlation-function scheme
In the relativistic quantum field theories, the Hilbert space only plays a secondary role, and we are supposed to reconstruct all the physical information from the correlation functions.
In the conventional Euclidean lattice gauge theory, people usually follow this dogma, and the mass spectrum is obtained from the correlation function in the imaginary time direction.We can take the same approach also in the Hamiltonian formalism by the measurement of the spatial correlation function with the distance r = |x − y|.
First, let us work on pions, and we consider the equal-time spatial correlation function, where π(x) denotes the operator defined by (3.24) with x = na.In order to evade the boundary effect as much as possible, we compute C π (r) by changing x and y symmetrically as x = (L − r)/2 and y = (L + r)/2, and the range of r is restricted to 0 ≤ r ≤ L/2.The results are shown in the left panel of Fig. 1 in the logarithmic scale, and the pion mass can be extracted from the exponential decay of C π (r).Here, the data with different colors represent the different values of the cutoff parameter ε.
It is convenient to use the so-called effective mass defined by where 2a comes from the step size of changing r.We further take the 3-point average of the effective mass One might be tempted to think that the pion mass corresponds to the plateau value of M π,eff (r), and the result with the DMRG cutoff ε = 10 −10 actually seems to become constant for r ≳ 10 almost exactly.However, this is the fake plateau due to the low-rank approximation.The point is that the leading asymptotic behavior of the spatial correlator is not purely the exponential decay, and it would take the Yukawa-type form asymptotically as r → ∞, We actually have α = 1/2 for the (1 + 1)d free massive boson, and we shall discuss the detailed analysis in Appendix D in the case of the 1-flavor Schwinger model.As a result, the effective mass for the Yukawa-type correlation function is given by and there must be an additional O(1/r) contribution on top of the actual mass M π .Motivated by this fact, we plot M π,eff (r) against 1/r in Fig. 2. We can see that the behavior of the effective mass strongly depends on the cutoff ε especially when r is large.When ε is not small enough, we observe the saturation of M π,eff (r), and then the α/r term seems to be absent.We note that the low-rank approximation of the DMRG is similar to the approximation of the transfer matrix by a finite matrix.Thus, C π (r) in the DMRG is approximated by the sum of purely exponential functions, and we need sufficiently large bond dimensions to reproduce the non-exponential corrections, such as 1/r α .
In fact, we can observe in Fig. 2 that the development of the 1/r-behavior in M π,eff (r) for large r by making ε sufficiently small, i.e. the bond dimension sufficiently large.We estimate the mass M π by the linear extrapolation 1/r → 0 of the result for the largest bond dimension with ε = 10 −16 , which is performed by fitting the data points with α/r + M π .To evaluate the systematic errors from the uncertainty of the fitting range, we try a lot of fittings by changing the fitting range inside the region 0.075 ≤ 1/r ≤ 0.15, and we obtain the probability distribution of the fitting results.The best-fitting result and its error are estimated from the position and the width of the peak, respectively.Thus, we obtained with α = 0.477 (9), and the fitting lines are drawn as the purple shadow in Fig. 2. Next, we perform similar analyses for sigma meson (3.26) and eta meson (3.25).Since these are isospin singlets, their one-point functions are not zero, and we subtract the disconnected parts from the correlation functions, with x = (L − r)/2, y = (L + r)/2.We then compute the 3-point averages of the effective mass, M σ,eff (r) and M η,eff (r), and they are shown in Fig. 3.The difference in the asymptotic behavior is observed by changing ε also in these cases.We plot the effective masses of the sigma and eta mesons against 1/r in Fig. 4 to see the asymptotic behavior.They approach ∝ 1/r for smaller ε as expected.We fit the data for ε = 10 −16 by α/r + M inside the region 0.075 ≤ 1/r ≤ 0.15 and estimate the systematic error.Then we obtained with α = 0.83(5) for sigma meson, and with α = 0.51(2) for eta meson.It is notable that α σ ∼ 0.8 has a relatively large deviation from the free boson result, α = 1/2, which may suggest that the sigma meson has a nontrivial dispersion relation even for small momentum.Finally, we summarize the masses of the three mesons measured by the correlation functions: pion sigma eta mass/g 0.431(1) 0.722(6) 0.899 (2) (5.11) The numerical results are qualitatively consistent with the analytic result by bosonization M π < M σ < M η .We also find M σ /M π = 1.68(2), (5.12) which is close to the prediction by the sine-Gordon model Eq.(2.19).

One-point-function scheme
We consider an alternative way to obtain the mass spectrum without using the two-point correlation functions.Let us recall that we are taking the open boundary condition, and we can use those boundaries as the source for excitations from the thermodynamic ground state.The boundary effect decays exponentially for the gapped systems, and thus the onepoint function of a local operator O(x) should behave as ⟨O⟩ + C ′ e −M O x as the function of the distance x = an from the boundary.Here, ⟨O⟩ gives the vacuum expectation value in the thermodynamic limit, and M O in the exponent gives the lightest particle mass with the same quantum number of O(x).In the context of condensed matter physics, it is known that the correlation function can be obtained from the Fridel oscillation, which is induced by a boundary effect or a local external field [4,5].We note that the x-dependence in this method takes the purely exponential form e −M x as the leading behavior for x → ∞.This can be easily understood by considering the path integral and the π/2 rotation of Euclidean spacetime.Then, the boundary condition sits at the constant imaginary time and defines the state |Bdry⟩ with zero momentum.Thus, the leading contribution to the imaginary-time correlation function ⟨Vac| Oe −H|x| |Bdry⟩ should come from the lightest particle with the zero-momentum projection, giving e −M x .This feature has nice compatibility with the low-rank approximation of DMRG.

5.2.1
The one-point functions of σ and η at θ = 0 At θ = 0, the boundary condition turns out to be completely invariant under the isospin rotation, and thus the boundary state |Bdry⟩ does not produce one pion states.Therefore, let us here focus on the iso-singlet particles, σ and η, and we will come back to pions later.
First, we discuss the eta meson as it turns out to be the simplest one.Since the G-parity is not spontaneously broken, we must have ⟨η⟩ = 0 in the thermodynamic limit.However, the staggered fermion realizes the G-parity (or charge conjugation) as the one-unit lattice translation, and thus the open boundary condition violates the G-parity.Therefore, the boundary state can be a source of the eta meson, and we evaluate the one-point function ⟨η(x)⟩ of the eta meson operator (3.25) in the range 0 < x ≤ L/2.The result is shown in with C = −1.096(1)for the smallest cutoff ε = 10 −16 .The errors of these values come from the fitting error.The corresponding fitting curve is shown in Fig. 5 with the purple line.In this case, we also find that the results for the other values of ε are consistent within the fitting error.Thus, the cutoff dependence does not appear unlike the case of the correlation functions, and we suspect that this is because MPS can efficiently express purely exponential decay.Next, we evaluate the one-point function ⟨σ(x)⟩ of the sigma meson (3.26) for 0 < x ≤ L/2.We note that σ has the same quantum number with the vacuum, and then ⟨σ(x)⟩ is nonzero also in the bulk.It behaves as e −M x+C + A with a constant shift of A, so we subtract the value ⟨σ(L/2)⟩ at the center x = L/2 of the lattice from ⟨σ(x)⟩ to remove the constant.The result is shown in Fig. 6, which indicates the exponential decay as expected.We fit the data points of ln | ⟨σ(x) − σ(L/2)⟩ | by −M σ x + C in the range 7 ≤ x ≤ 13, and the best-fit parameter is M σ = 0.761(2), (5.14) with C = −2.71(2), which are independent of the value of ε up to the fitting error.The result of fitting for ε = 10 −16 is shown in Fig. 6 with the purple line.

The one-point functions of π 3 at θ = 2π
Let us now come back to the issue of pions.As we have argued, the boundary state at θ = 0 is neutral under the isospin rotation, and thus it does not produce one-pion states and we have ⟨π(x)⟩ = 0 for all x.Therefore, we need to somehow create the boundary state that transforms nontrivially under the isospin rotation to study pions with the onepoint-function scheme.In this study, we decided to use one of the ground states at θ = 2π for this purpose.Since the Hamiltonians at θ = 0 and θ = 2π are unitary equivalent under the periodic boundary condition, the bulk properties are the exactly same between θ = 0, 2π.As we have discussed in Section 2.2, the ground state at θ = 2π is a nontrivial SPT state protected by the isospin SU(2) V /Z 2 symmetry, and thus the boundary states with the open boundary condition have the isospin 1/2.This boundary charge can be a source of the pions so that ⟨π(x)⟩ becomes nonzero.
About the computational cost, it turns out that the bond dimensions for the MPS are mostly the same at θ = 0 and θ = 2π when the system size is large enough.Therefore, we can obtain the ground state at θ = 2π as easily as that of θ = 0. We, however, observe that the bond dimension at θ = 2π increases significantly if the system size is not large enough, and we suspect its reason is as follows.At θ = 2π, there is 4-fold degeneracy due to the boundary degrees of freedom, but they split into the singlet and the triplet states with the energy splitting ∼ e −MπL .That is, the true ground state at finite L has an extra Bell pair between the boundary isospin 1/2 states, which adds ln 2 to the entanglement entropy.When we cut the system at x = L/2, this extra ln 2 should be accumulated by the large numbers of small singular values, and thus the bond dimension becomes quite huge just to create the Bell pair between the boundaries.If L is large enough, the energy gain by creating the Bell pair becomes negligible, and the DMRG would produce one of the ground states with disentangled boundary states practically.Thus, the computational cost becomes almost the same as that for the trivial state at θ = 0.
Let us now evaluate the one-point function ⟨π(x)⟩ of the pion (3.24) for 0 < x ≤ L/2 using the ground state at θ = 2π.The result is shown in Fig. 7.We again find the exponential decay, and thus fit the data points of ln with C = 0.203 (9), which do not depend on the cutoff ε up to the fitting error as before.The fitting result for ε = 10 −16 is shown in Fig. 7 with the purple line.
Let us summarize the effective masses using the one-point-function scheme: pion sigma eta mass/g 0.4175(9) 0.761(2) 0.9014(1) ( The order of three meson masses is consistent with the analytic prediction.We also find which is still close to the WKB prediction, √ 3, with a 5% deviation.The significant feature of the method is that the results do not depend on the cutoff parameter ε as long as it is sufficiently small.Therefore, the systematic error from the cutoff is expected to be small enough.We do not need to increase the bond dimension so much, unlike the method by the correlation function.

Dispersion-relation scheme
So far we have studied the mass spectrum by using the local observables of the ground state, and these methods are applicable both in the path integral and the Hamiltonian formalisms.As the third method for computing the mass spectrum, let us take a different approach that is specific to the Hamiltonian formalism: We compute the excited states as explained in Section 4.1, and then determine the mass spectrum from the dispersion relation.
The low-lying excited states correspond to one-particle excitations.For example, the zero mode of the lightest meson, namely the pion, is expected to be obtained as the first excited state.We can also obtain the states with nonzero momentum K, and we can fit the data with the dispersion relation ∆E ≃ K 2 + M 2 π to obtain the pion mass.As we go to the higher excited states, we will encounter one-particle states of the sigma and eta mesons.They can be distinguished by measuring quantum numbers, such as the isospin and G-parity.Thus, we can compute the mass spectrum from the dispersion relation by generating the excited states.
We note that our computation is done in the finite open interval, and thus the momentum is not a good quantum number.Also, there may exist a nontrivial contribution to the excitation energy from the boundaries.We are neglecting those subtleties in this work, but, surprisingly, it turns out that the numerical results are almost consistent with those with the previous two methods.
We generated the MPS up to the 23rd excited state at θ = 0 with the small physical volume L = 19.8.The energy gap ∆E ℓ = E ℓ − E 0 of the ℓ-th excited state is shown in the left panel of Fig. 8.We also measured the square of total momentum K 2 defined by (3.46).We note that its ground-state expectation value K 2 0 ≃ 0.46 is nonzero because of the boundary effect and maybe also due to lattice artifacts, and thus we subtract K 2 0 from K 2 ℓ of the excited states.The result is plotted in the right panel of Fig. 8. From these results, we find many triply degenerated states, which are candidates for the states of the pion.There are a few singlet states as well, which are candidates for the eta and sigma mesons.
To identify the states, we measure the expectation values of the isospin operators, J 2 and J z , the parity P and the G-parity G = Ce iπJy defined in Section 3.4.We note that the DMRG does not produce the states in a diagonal basis for these quantities.We diagonalize the 3 × 3 matrix ⟨ψ ℓ 1 | J z | ψ ℓ 2 ⟩ in each triplet to compute expectation values in the J z basis. 6,7 The expectation values of J 2 , J z , G, and P in the J z basis are listed in Tables 1  and 2 for iso-triplets and iso-singlets, respectively.The index ℓ comes from the level of the state on the original random basis.We find that |G| ̸ = 1 because of |C| ̸ = 1 by the effect of the boundary.Hopefully, the sign of G can be assumed to remember the original ℓ J 2  J z G P  quantum number [25], and, if it is true, we can still identify the G-parity.This point will be discussed more in details in Appendix B. We identify the lowest triplet ℓ = 1, 2, 3 as the lowest modes of the pions (π + , π 0 , π − ) since they have the quantum numbers consistent with the pion, namely J P G = 1 −+ and J z = 0, ±1.For the iso-singlets shown in Table 2, we find that the ℓ = 13 state has the quantum number consistent with the sigma meson, namely J P G = 0 ++ and J z = 0.The ℓ = 18 state is consistent with the eta meson with J P G = 0 −− and J z = 0. We identify these singlets with the lowest modes of the sigma and eta mesons.
After identifying the quantum numbers, we plot the energy gap ∆E ℓ = E ℓ −E 0 against the momentum square ∆K  We find the mass ratio M σ /M π ≃ 1.75(1) (5.19) from this result, which is close to the WKB prediction √ 3.

Conclusion and Discussion
In this paper, we work on three independent methods to compute the mass spectrum of lattice gauge theories in the Hamiltonian formalism, which apply to tensor networks and quantum computation.The methods are tested in the massive 2-flavor Schwinger model at θ = 0, some of which properties are analogous to the ones of 4d QCD.The two species of fermion play roles of up and down quarks, and the composite particles (mesons) appear as triplets or singlets of the SU(2) V /Z 2 isospin symmetry.We used the tensor network, in particular, DMRG for numerical simulation. .The masses of the pion, sigma, and eta meson obtained by the three independent methods are compared.Each result is obtained with the given finite lattice spacing.We also put the error bar of the fitting error for the correlation-function scheme, but it is too small to be seen.
We obtained the masses of the pion, sigma, and eta meson by the three methods, and the results are summarized in Fig. 10.We find that the results are roughly consistent with each other taking into account possible systematic errors for each method, such as the continuum and infinite-volume limits.In addition, all the results show the relation M π < M σ < M η , which agrees with the analytic prediction by the bosonization technique.The order of the eta meson mass M η ∼ 0.9 is consistent with M η ∼ µ since µ ∼ 0.8 in the current setup.We also find that the relation between the masses of the pion and sigma mesons is M σ /M π = 1.68(2), 1.821 (6), 1.75(1) by the correlation-function scheme, the one-point-function scheme, and the dispersion-relation scheme, respectively.We note that the errors in the above values only contain the fitting error, and there should be further systematic errors potentially coming from the finite lattice spacing, the finite-volume effect, the open boundary condition, the cutoff of the bond dimension, etc.These results are close to the WKB-based formula (2.19), M σ /M π = √ 3, within not more than a 5% deviation.It is, honestly, very surprising that the semiclassical analysis gives the almost correct answer outside the range of its validity, and it would be theoretically interesting to uncover the reason behind its success.
Let us discuss the advantages and difficulties of each method and the potential applications to other models.The first one, the correlation-function scheme, is the straightforward generalization of the technique in Lagrangian formalism.The advantage of this method is a wide range of applicability to various models.We can obtain the meson masses from correlation functions on a lattice with any dimensions, volume, and boundary condition.Furthermore, the correlation function accepts the off-diagonal element such as ⟨O(x)O ′ (y)⟩.This feature will be useful when we turn on θ ̸ = 0 in the 2-flavor Schwinger model.The reason is that the meson operators become nontrivial mixtures of S f (n) and P S f (n) depending on θ.In this case, we need to measure the correlation matrix of the operators and diagonalize it to extract the mode of each meson.However, our numerical results suggest that the bond dimension of MPS has to be sufficiently large to reproduce the correct asymptotic behavior of the correlation function.In particular, the computational cost increases rapidly as the system approaches a gapless phase, for example, m ∼ 0 or θ ∼ π.
Thus, the tensor network (MPS) is not an efficient approach to computing the mass spectrum by using correlation functions. 9On the other hand, an ideal quantum computer is free from such a restriction of the bound dimension.Thus, the correlation function may be the first option in the era of practical quantum computation of field theories in this sense, though to avoid the finite volume effect for the two-point function we need a sizable scale computer.
The second method, the one-point-function scheme, makes good use of the boundary effect rather than eliminating it.The results turn out to be insensitive to the bond dimension, and thus we have to increase neither the lattice size nor the bond dimension so much.Furthermore, the evaluation of the local one-point function is generally easier than that of the long-range correlation function.Thus, this is the most economical one among the three methods.We note, however, that we have to prepare suitable boundary conditions, such as defects, impurities, or external fields, to compute the mass spectrum with this one-point-function scheme, which requires good physical insights for the system of interest.In our case, the open boundary at θ = 0 can be regarded as a source of the iso-singlet mesons, σ, and η, but we have to set θ = 2π to induce the boundary excitation as a source of the iso-triplet mesons, π a .We should also note that we cannot obtain information on the off-diagonal correlators in the one-point-function scheme.When θ = 0, the off-diagonal correlators are unimportant because π, σ, and η have different quantum numbers, but they should become important at generic values of θ because the G-parity is no longer a good quantum number.
The third method, the dispersion-relation scheme, is the distinctive strategy of Hamiltonian formalism.We can obtain various states heuristically without knowing what kind of mesons appear in the spectrum.Once we generate the excited states, it is straightforward to measure various observables such as energy, momentum, and quantum numbers.The states are identified by using these pieces of information.Furthermore, we can investigate the wave function to distinguish the s-wave or p-wave states.In this method, however, it is difficult to increase the system size or the spatial dimensions.The reason is that we have to generate an increasing number of states to search for different mesons.For example, in our setup, we encounter the 3 × 4 states of the pion before obtaining the sigma meson at ℓ = 13.The momentum K is discretized as K ∼ 2πκ/L for κ = 1, 2, • • • in the finite system with the size L. If L is increased, the number of pion states in a certain range of energy grows up.Thus, we have to generate more excited states to reach the state of the sigma meson.As for higher dimensions, there are momentum excitations in each spatial direction, which result in an additional degeneracy.We expect that there is a way to avoid this issue by modifying the strategy.For example, if we are interested in a specific meson, it is more effective to generate excited states with a constraint on the quantum number to skip mesons out of interest.
In this work, we have computed the mass spectrum at θ = 0. We note that we have neglected many systematic errors, and thus there is plenty of room for improvement.As a physics, extending our investigation to θ ̸ = 0 should be interesting, where the sign problem arises in naive applications of Monte Carlo simulations.The presence of θ introduces some differences compared to the θ = 0 case.Firstly, the mass of the pion, which corresponds to the gap of the system, decreases as θ → π.Consequently, we may need to increase the bound dimension of MPS, leading to higher computational costs.Secondly, the parity and G-parity are no longer good quantum numbers for θ ̸ = 0, and the scalar and pseudoscalar operators have a nontrivial mixture.To handle this situation, we should measure the correlation matrix between these operators and diagonalize it.Although distinguishing the excited states, especially σ and η, seems to become tricky, exploring the changes in the spectrum promises intriguing insights.Despite these subtleties, we expect that it is still possible to apply the three methods to compute the mass spectrum including the theta term, and the results at θ ̸ = 0 will be reported elsewhere.Needless to say, it is very desirable that future developments of these techniques eventually enable us to compute the hadron spectrum of 4d strongly-coupled gauge theories having the sign problem in the conventional Monte Carlo methods.First, we investigate the behavior of C in the continuum limit.We generated the MPS of the ground state and the 1st excited state of the 1-flavor Schwinger model at θ = 0 by DMRG.The lattice spacing a is changed around 0.1 ≲ a ≲ 0.25.The number of lattice sites N is chosen to fix the physical system size L = (N −1)a.We compute the expectation values of C for these MPS.The results are shown in Fig. 11.The different symbols correspond to the results for different L in the plots.We fitted the data points for each L by the quadratic function f (a) = c 0 + c 1 a + c 2 a 2 .The fitting results are also plotted in Fig. 11 by the solid lines.For L = 49.8,we obtained the continuum limit ⟨C⟩ a→0 = 0.321(3) for the ground state and −0.320(3) for the 1st excited state.The results with the other L agree with these values within the error.Thus, we confirmed that the expectation value of C is a finite value in the continuum limit and is not sensitive to L.
Next, let us discuss the effect of the boundary on C. We consider a further simplified model, the free fermion on the periodic lattice.The model is obtained from the 1-flavor Schwinger model with the periodic boundary condition by setting g = 0 and adding the hopping term between n = 0 and n = N − 1 site.In fact, it is hard to adopt the p.b.c. in the current DMRG method due to the artificial long-range interaction between both ends of MPS.Thus, we choose small sizes of the lattice N = 20 and 40 for this analysis.The corresponding lattice spacings are a = 0.2 and 0.1 for the fixed physical length L = N a = 4.We generate the ground state and the excited states up to the level ℓ = 4.The four excited states turned out to be degenerated.Thus, we compute ⟨C⟩ ℓ,ℓ ′ including the off-diagonal elements, and diagonalize the result as the 4 × 4 matrix.The eigenvalues are shown in Fig. 12.We found that ⟨C⟩ = 1 for the ground state and ⟨C⟩ = ±1, α ± iβ for the excited states.These complex values satisfy |α| 2 +|β| 2 = 1 as we can see in the plot.The imaginary part β becomes smaller as a is decreased, which suggests that we will obtain ⟨C⟩ → ±1 in the continuum limit.which we name the flavor order here.In this case, we first arrange one of the flavors and then start to arrange the next one, so the flavor degrees of freedom at the same physical sites are separated by N , and this clearly violates the above important criterion.In these two cases, we compare the efficiency of the MPS to represent the ground state in the gapped phase θ = 0. We obtain the ground state by DMRG and investigate the largest bond dimension in the MPS, called the effective bond dimension D eff .The results are plotted against the number of sweeps N sweep in Fig. 13 for various lattice sizes N .We found that D eff converges around O(10) sweeps for both cases.However, the dependence on N is totally different.For the flavor order, the final value of D eff increases exponentially with N , which is caused by artificial long-range interaction between the two flavors.On the other hand for the staggered order, the final value is saturated for sufficiently large N .To show these behaviors, we plot the final values of D eff after 20 sweeps against N in Fig. 14.According to Fig. 14, the bond dimension seems to saturate in the case of the staggered order as N → ∞.Since the ln D eff gives the upper bound for the entanglement entropy, this constant behavior is expected to be the optimal one for the 1 + 1d gapped systems.On the other hand, D eff grows exponentially fast for the flavor ordered as N → ∞.We suspect that this is because the flavor order puts the entangled flavors in separate locations.If we cut the system into two pieces in terms of i with the flavor ordering, the O(N ) entangled pairs are cut, and thus the entanglement entropy becomes O(N ), which is consistent with the exponential behavior of D eff .Therefore, we adopt the staggered order in the whole analysis of this work.

D Correlation function in the 1-flavor Schwinger model
We test the validity of the correlation-function scheme in Section 5.1 by examining the correlation function in the 1-flavor Schwinger model.When the fermion is massless m = 0, the model can be analytically solvable and it is equivalent to the free massive boson with mass µ ′ = g/ √ π.Thus, this is a good benchmark and we compare the numerical result of DMRG with the analytical answer.
As an analogy of the pseudo scalar meson in the 2-flavor Schwinger model, we consider the pseudo-scalar operator P S = −i ψγ 5 ψ.The results of the correlation function ⟨P S(x)P S(y)⟩ are shown in the left panel of Fig. 15.Here, the data with different colors are obtained with the different values of the cutoff parameter ε.The corresponding effective masses (3-point average) are also plotted in the right panel of Fig. 15, where we can see the significant ε dependence.
To see the 1/r correction of the effective mass, we plot M eff (r) against 1/r in Fig. 16.Then we found that the result approaches the expected asymptotic behavior M eff (r) ∼ α/r + M only if the cutoff ε is sufficiently small.We fitted the data points for ε = 10 −16 by α/r+M in the range 0.06 ≤ 1/r ≤ 0.2 and obtained M = 0.5677 (5) and α = 0.446 (4).Here the systematic error from the uncertainty of the fitting range is evaluated as explained in Section 5.1.We note that this is the result on the finite lattice before taking the continuum limit, but it turned out to be close to the exact value M = g/ √ π ≈ 0.56419 of the continuum theory.
Therefore, it is quite important to discuss the cutoff (or bond-dimension) dependence especially when we use the correlation-function scheme.Indeed, if we naively read the plateau value at ε = 10 −10 , we got an incorrect answer M ∼ 0.63 without observing the 1/ √ r contribution in the Yukawa potential at all.

A
Operators in the spin representation B Charge conjugation operator in the 1-flavor Schwinger model C Arrangement of flavors on MPS D Correlation function in the 1-flavor Schwinger model 1 Introduction

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

Figure 1 .
Figure 1.(Left) The correlation function of pion ln | ⟨π(x)π(y)⟩ | is plotted against the distance r = |x − y| for various values of ε.The number of lattice sites is N = 160 and the lattice spacing a is determined so that L = a(N − 1) = 39.8.(Right) The effective mass of the pion M π,eff (r) (3-point average) calculated from the correlation function in the left panel is plotted against r.

MFigure 2 .
Figure 2. The effective mass of the pion M π,eff (r) is plotted against 1/r.The data points for ε = 10 −16 are fitted by α/r + M inside the region 0.075 ≤ 1/r ≤ 0.15.The fitting result is depicted by the shaded band with systematic error.

Figure 4 .
Figure 4.The effective mass of sigma meson M σ,eff (r) (left) and eta meson M η,eff (r) (right) are plotted against 1/r.The data points for ε = 10 −16 are fitted by α/r + M inside the range 0.075 ≤ 1/r ≤ 0.15.The fitting results are depicted by the shaded bands with systematic errors.

Figure 5 .
Figure 5.The one-point function ln | ⟨η(x)⟩ | of the eta meson is plotted against x = an with n = 1, • • • , N/2 − 1 for various values of ε.The number of lattice sites is N = 160 and the lattice spacing a is determined so that L = a(N − 1) = 39.8.The result of fitting by −M η x + C for ε = 10 −16 is also plotted by the solid line inside the range and by the broken line outside.

Fig. 5 .
Fig. 5.The cutoff parameter is changed from ε = 10 −10 to 10 −16 .The one-point function decays exponentially with x as expected.Thus, we fit the data points of ln | ⟨η(x)⟩ | by −M η x + C in the fitting range 7 ≤ x ≤ 13, and the result is

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

Figure 7 .
Figure 7.The one-point function ln | ⟨π(x)⟩ | of the pion is plotted against x = an with n = 1, • • • , N/2 − 1 for various values of ε.The number of lattice sites is N = 160 and the lattice spacing a is determined so that L = a(N − 1) = 39.8.We set θ = 2π in order to induce the boundary charges, which make ⟨π(x)⟩ nonzero.The result of fitting by −M π x + C is also plotted by the solid line inside the range and by the broken line outside.

Figure 8 .
Figure 8. (Left) The energy gap ∆E ℓ = E ℓ − E 0 is plotted against the level of the excited state ℓ.(Right) The square of total momentum ∆K 2 ℓ = K 2 ℓ − K 2 0 is plotted against ℓ after subtracting the result for the ground state.

26 -Figure 9 .
Figure 9.The energy gap ∆E ℓ is plotted against the square of total momentum ∆K 2 ℓ .The states with the same isospin and G-parity are plotted by the same symbol.Then each state is identified with the pion, sigma, or eta meson.We fit the data for each meson by ∆E = √ b 2 ∆K 2 + M 2 .The results are shown by the broken lines.The values of M for each meson are also plotted as the endpoints of the fitting lines.

Fig. 9 .
Fig.9.The states with the same isospin J 2 and G-parity are plotted by the same symbol.8Then we fit the data points by ∆E = √ b 2 ∆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 ∆K 2 → 0. We obtained M π = 0.426(2), b π = 1.017(4) for the pion; and M σ = 0.7456(5), b σ = 1.087(2) for the sigma meson with the fitting error.The fitting for the eta meson is simply solving an equation since there are only two data points.The result are M η = 0.904 and b η = 0.962.We summarize the masses of the mesons determined by the energy gap of the excited states:

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

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

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

Figure 13 .DFigure 14 .
Figure 13.The effective bond dimension D eff is plotted against the number of sweeps N sweep for the flavor order (left) and the staggered order (right).The vertical axis of the left panel is in log scale, whereas the axis of the right panel is in linear scale.The lattice spacing and the fermion mass are set to a = 0.2 and m = 0.1.

Figure 15 .
Figure 15.(Left) The correlation function ln ⟨P S(x)P S(y)⟩ is plotted against the distance r = |x − y| for various values of ε after subtracting the disconnected part.The number of lattice sites is N = 400 and the lattice spacing a is determined so that L = a(N − 1) = 79.8.(Right) The effective mass M eff (r) (3-point average) calculated from the correlation function in the left panel is plotted against r.

Figure 16 .
Figure16.The effective mass M eff (r) is plotted against 1/r.The data points for ε = 10 −16 are fitted by α/r + M inside the region 0.06 ≤ 1/r ≤ 0.2.The fitting result is depicted by the shaded band with systematic error.The exact mass of the pseudo scalar g/ √ π is also shown by the horizontal broken line.

Table 1 .
The quantum numbers of the isospin triplet states.The index ℓ comes from the level of each state in the original basis.The rows of the table are separated into each triplet.

Table 2 .
The quantum numbers of the isospin singlet states.