The mass spectrum of the Schwinger model with matrix product states

We show the feasibility of tensor network solutions for lattice gauge theories in Hamiltonian formulation by applying matrix product states algorithms to the Schwinger model with zero and non-vanishing fermion mass. We introduce new techniques to compute excitations in a system with open boundary conditions, and to identify the states corresponding to low momentum and different quantum numbers in the continuum. For the ground state and both the vector and scalar mass gaps in the massive case, the MPS technique attains precisions comparable to the best results available from other techniques.


JHEP11(2013)158
A particularly interesting study case is that of the Schwinger model [25,26], or QED in one spatial dimension, the simplest gauge theory, which nevertheless exhibits features in common with more complex models (QCD) such as confinement or a non-trivial vacuum, and has been adopted as a benchmark model where to explore lattice gauge theory techniques (see e.g. [27][28][29][30][31] and references cited therein).
The application of MPS techniques to the Schwinger model was first explored by Byrnes et al. [20] using the original DMRG formulation. The study improved by several orders of magnitude the precision attained by other Hamiltonian techniques for the ground state and the first (vector) mass gap, for vanishing and non-vanishing fermion masses, using periodic boundary conditions (PBC). However, DMRG estimations for higher excited states lose precision fast, in particular for PBC, and the scalar mass gap was not explored in [20].
In this paper we apply the more stable and numerically efficient MPS for open boundary conditions (OBC) to the same problem. We work in the subspace of physical states satisfying Gauss' law, in which the model can be written as a spin Hamiltonian with a long-range interaction [32,33]. We improve the existing techniques to find higher excited states and devise a method to identify vector and scalar excitations in finite open chains, although the charge that distinguishes them is only a good quantum number in the continuum. Using these methods we compute the ground state and the vector and scalar mass gaps for vanishing and non-vanishing fermion masses, with enough precision to conduct the extrapolation to the continuum limit. Our study shows the feasibility of TNS solutions for similar LGT, in the Hamiltonian formulation.
The rest of the paper is organized as follows. In section 2 we briefly introduce the Schwinger model, and review its formulation as a spin Hamiltonian. Section 3 presents the MPS formalism, with special emphasis on the new techniques used to overcome the challenge posed by this particular kind of problems. In section 4 we present our numerical results and conclude in section 5 with a discussion and outlook.

The model
The massive Schwinger model [25], or QED in two space-time dimensions, is a gauge theory describing the interaction of a single fermionic species with the photon field, via the Lagrangian density, with F µν = ∂ µ A ν − ∂ ν A µ , being g the coupling constant and m the fermion mass. The physics of the model is determined by the only dimensionless parameter m/g. The massless case, m = 0, can be solved analytically, via bosonization [25], and also the free case, g = 0, has exact solution, so that for very large or very small values of m/g a perturbative study is possible around one of the solvable limits, while in general a non-perturbative treatment may be needed. One of the features of the model is the existence of bound states, the two lowest ones of which, the vector and the scalar, are stable for any value of m/g [26,34]. The spectrum of the Schwinger model has been studied with many techniques. The most accurate numerical estimations have been performed on the lattice. For the vector JHEP11(2013)158 mass the best results were obtained with DMRG [20], while strong coupling expansion (SCE) got the most accurate prediction for the scalar mass [35,36] and the most precise values in the massless case [37].
In the temporal gauge, A 0 = 0, the Hamiltonian density reads The electric field, E = −Ȧ 1 , is fixed by the additional constraint of Gauss' law, ∂ 1 E = gΨγ 0 Ψ, up to a constant of integration, which can be interpreted as a background field [26]. The model can be formulated on a lattice. Here we focus on the Kogut-Susskind staggered formulation [38], where a is the lattice spacing. We consider a lattice with open boundaries, with sites n = 0, . . . N − 1. On each site there is a single-component fermion field φ n , while θ n are the gauge variables sitting on the links between n and n + 1, and connected to the vector potential via θ n = −aqA 1 n . L n is the corresponding conjugate variable, [θ n , L m ] = iδ nm , connected to the electric field as gL n = E n . We work in a compact formulation, where θ n becomes an angular variable, and L n can adopt integer eigenvalues. Gauss' law appears as the additional constraint Instead of the Hamiltonian (2.3), it is usual to work with the dimensionless operator W = 2 ag 2 H, and to define parameters x = 1 g 2 a 2 and µ = 2m g 2 a . In this picture, fermions on even and odd sites respectively correspond to the upper and lower components of the spinor representing the fermionic field in the continuum. Gauss' law (2.4) determines the electric field up to a constant α which can be added to L n and represents the background electric field.

JHEP11(2013)158
where ℓ is the boundary electric field, on the link to the left of site 0, which can describe the background field. A useful basis for this problem is then with i m = {0, 1} labeling the ±1 eigenstates of σ z m (on site m). An even site in spin state |1 corresponds to the presence of a fermion, while an odd site in state |0 is an antifermion at the corresponding position. The integer valued ℓ = . . . , −1, 0, 1, . . . is the only gauge degree of freedom left, but with OBC it is non-dynamical, as the Hamiltonian cannot connect states with different values of ℓ. Here we choose ℓ = 0 and omit it in the following, as we will be concerned with the case of zero background field. Moreover, we are interested in states with zero total charge ( n σ z n = 0 in the spin language), so we consider chains with even N .

Method
In this work we use the MPS ansatz to approximate the ground and lowest excited states of the Hamiltonian (2.6). A MPS for a system of N d-dimensional sites has the form are individual basis states for each site. Each A i k is a D-dimensional matrix and the bond dimension, D, determines the number of parameters in the ansatz. MPS are known to provide good approximations to ground states of local Hamiltonians in the gapped phase, but have also been successfully used for more general models.
The MPS approximation to the ground state can be found variationally by successively minimizing the energy, Ψ|H|Ψ Ψ|Ψ , with respect to each tensor A k until convergence. Each such optimization reduces to the eigenvalue problem of an effective Hamiltonian that acts on site k and its adjacent virtual bonds. The basic algorithm, closely related to the original DMRG formulation [7,8], was introduced in [5], and is nowadays widely used.
In DMRG it is possible to target several of the lowest eigenstates [41], or, in the case of a first excited state with a distinct quantum number, to run the ground state search in different sectors. In the MPS formalism, it is natural to extend the variational method to determine excited states by performing a constrained minimization which imposes orthogonality with respect to the already computed states (ground and lower excited states), as proposed in [42] and recently extended in [43]. Here we introduce a simpler variation of the ground state algorithm, which allows us to compute MPS approximations to the lowest energy eigenstates at lower cost, without using any explicit symmetry. This is especially interesting for a finite system with open boundary conditions, where momentum is not a good quantum number. It is worth noticing that in the case of translationally invariant (TI) systems, either finite and periodic or infinite, a very successful approach exists based on the construction of well-defined momentum MPS [44].

JHEP11(2013)158
For finite systems, the MPS algorithms for open boundary conditions are numerically more stable and in general more efficient, although improved methods have been recently proposed for periodic systems [45,46]. However, finite-size effects are much larger for OBC and therefore simulation of larger chains may be needed to reach the thermodynamic limit reliably.
After having found the ground state of the system, |Ψ 0 , we can construct the projector onto the orthogonal subspace, Π 0 = 1 − |Ψ 0 Ψ 0 |. The projected Hamiltonian, Π 0 HΠ 0 , has |Ψ 0 as eigenstate with zero eigenvalue, and the first excited state as eigenstate with energy E 1 . Given that E 1 < 0, what we can always ensure by adding an appropriate constant to H, the first excitation corresponds then to the state that minimizes the energy of the projected Hamiltonian, 1 This minimization corresponds to finding the ground state of the effective Hamiltonian H eff [1] = Π 0 HΠ 0 . The procedure can be concatenated to find subsequent energy levels, so that, to find the M -th excited state, we will search for the ground state of the Hamiltonian Each of these ground state searches can be solved by applying the standard variational MPS algorithm to the corresponding effective Hamiltonian (3.3), which can be constructed from the set of all previously computed MPS levels {|Ψ k } (see figure 2). The expression to be minimized at each step of the MPS iteration is then min A k where H k and N k are the effective Hamiltonian and norm matrix acting on site k and its two adjacent virtual bonds, so that each tensor can be found by solving a generalized eigenvalue problem H k A k = λN k A k . As illustrated in figure 1, H k and N k are computed by contracting all tensors but A k , and the computational cost can be optimized, as in the original algorithm, by storing and reusing the partial contractions that compose these effective matrices. In general this is most easily done when all the terms in the problem Hamiltonian are expressed as a matrix product operator (MPO) [47]. In this particular case it is more convenient to keep separate temporary terms for each of the contractions Ψ k |Ψ , and to reconstruct from them the effective projectors on every site, before constructing the effective Hamiltonian. In this way, the number of operations required for each such term scales as dD 3 , without increasing the leading cost of the algorithm. Due to the larger number of terms that need to be kept, however, the cost of finding the M -th level once all the lower ones are known will be about M times higher than that of the original ground state search, 2 and thus the total cost for computing up to the M -th level will scale as M 2 D 3 .
The masses of the particles in the theory are given by the energies of the zero momentum excitations. In finite systems with OBC momentum is difficult to identify and we   need to find the excitations that will correspond to the lowest momentum in the continuum limit. In dynamical DMRG [48] momentum dependent quantities were extracted from finite DMRG calculations with open boundaries using quasimomentum states defined from the eigenstates of a free particle in a box. Here we use a different approach, based on the continuum momentum operator for the fermion field,P = dxΨ † (x)i∂ x Ψ(x). Its discretization yields, in the spin representation, and after rescaling by a factor 2 ag 2 , the operator The expectation values ofÔ 2 P can be used to assign a pseudo-momentum to the spectral levels, in order to reconstruct the dispersion relation (figure 3).

JHEP11(2013)158
Input: Hamiltonian for an N -site chain, H; maximum bond dimension, D; tolerance, ǫ Output: MPS approximation to the ground state,   The appearance of the scalar is detected by the phase of the expectation value S R , with the scalar (indicated in blue in the plots) being the first state with ϕ ≈ 0 above the ground state (red) while states with |ϕ| ≈ π belong to the vector branch (green).

JHEP11(2013)158
The model gives rise to stable particles, the lowest ones being, in the case of no background field, the vector and scalar states. In the continuum model, these particles are distinguished by well-defined parity and charge conjugation quantum numbers [26], with the scalar living in the same sector as the ground state. On a staggered lattice with PBC it is possible to exploit the corresponding lattice symmetries to construct two orthogonal subspaces, one of them containing the vector, and the other the ground state and the scalar [49]. For a chain with OBC the number of surviving symmetries is even lower, with translational invariance lost. In practice this means that to calculate the scalar mass, we need to identify first the momentum excitations of the vector, which appear at lower energy than the scalar. In [39] this was achieved by starting from the strong coupling limit (x → 0), where vector and scalar states are known exactly, and smoothly changing the x parameter while keeping a label on the scalar state. Instead, we use the expectation value of the spin transformation, is the (cyclic) translation by one spin site. 3 In a system with PBC this operator basically describes the action of charge conjugation on the spins, with the translation exchanging the fermionic and antifermionic character of sites, and the σ x rotation accounting for the different meaning of a spin up on an even (fermion) site (empty) and on an odd one (occupied). Charge conjugation is a good quantum number and distinguishes the vector state (C = −1) from the sector containing the scalar and the ground state (C = +1), even for finite systems. On the staggered lattice with OBC this is no longer true, and S R does not commute with the Hamiltonian. However the phase of S R keeps memory of this distinction, and allows us to tag the levels accordingly, with the ground state and the scalar branch of the dispersion relation having phase ϕ( S R ) ≈ 0 and the vector branch ϕ( S R ) ≈ π, as illustrated in figure 3.

Results
We have applied the methods described in section 3 to the Hamiltonian (2.6) to determine the ground state, and the vector and scalar mass gaps for various values of the fermion mass. In order to benchmark our results with existing data, we have studied different masses m/g = 0, 0.125, 0.25 and 0.5, for which reference values can be found in the literature.
Our goal is to compute, for each of these masses, the ground state energy density, and the vector and scalar mass gaps in the thermodynamic limit (N → ∞) for different finite values of x, and to extrapolate them to the continuum limit, corresponding to the lattice spacing a → 0, or x → ∞. We have used values of x ∈ [5,600]. In order to extract the thermodynamic limit for each case, we need to simulate different system sizes, N , and apply finite-size scaling. The system sizes, N , cannot be chosen independently of the value of x, as the same number of sites N corresponds to different physical volumes L phys = N a for different values of x ∝ 1/a 2 . From numerical simulations, we estimated that N ≥ 20 √ x ensures small enough finite volume effects. We thus study for each value of x several different system sizes, large enough for the condition above to be satisfied (for instance, for x = 5 we take N ∈ [50,82] while for x = 600, N ∈ [540, 834]).

JHEP11(2013)158
We have run the MPS algorithms described in the previous section for each set of parameters, (m/g, x, N ), using bond dimensions D ∈ [20,140], to find the ground state and excited levels at least until a candidate scalar state is found. We are interested in the subspace of null total charge, which corresponds to S z = 0 and can easily be imposed by adding a penalty term to the Hamiltonian. For every level, the MPS iteration stops when the relative change in the energy after one full sweep over the chain is below a certain tolerance, ǫ. This value has to be small to ensure a good precision after the extrapolations, but a smaller tolerance translates in more difficult convergence. Therefore we fixed ǫ = 10 −12 for the ground state and the vector calculations, and ǫ = 10 −7 for the longer computations required by the scalar states.
For every set (m/g, x, N ) and for each of the levels we are interested in, we extrapolate our results to 1/D → 0, as illustrated in figure 4. In the case of the ground state and the vector state, almost all the bond dimensions are converged, while for the scalar mass gap, depending on a larger number of intermediate excited states, some of the larger systems are only computed with smaller bond dimension. If the bond dimension reached is large enough, we extrapolate linearly in 1/D. Otherwise we take the value corresponding to the largest D as our estimation for the energy, and the error as the difference between this value and the one for the immediately smaller D.
The results of the extrapolation described above provide accurate estimators for the energy levels for various finite chains. We then proceed to scale these results with the finite system size to extract the ground state energy density, E 0 /(2xN ), and the mass gaps, (E 1(2) − E 0 )/(2 √ x) in the thermodynamic limit for each pair of parameters (m/g, x).
Finite-size corrections to the ground state energy density are linear in 1/N , while for the energy gaps, the leading corrections arise from a kinetic energy term O(π/N 2 ) [39]. Hence the bulk limit is extracted by fitting the ground state energy density as and the energy gaps as Finally, the values in the continuum limit can be extracted by fitting the results for each value of the mass to a polynomial in 1 √ x . We include only those values of x for which the thermodynamic limit could be extracted accurately (i.e. the corresponding level for the various system sizes was computed for a large enough D). The interval [x min , x max ] and the degree of the polynomial used for the fit introduce the largest uncertainty in the final estimators. We have thus performed a systematic error analysis in which the different possible fits are taken into account to estimate the continuum values and their errors from our data. The detailed method is described in the appendix A.  As our final result we obtain for the ground state energy density in the massless case the value ω 0 = −0.318338 (22) (24), to be compared to the exact value, −1/π ≈ −0.3183099. The first error includes the propagated error from infinite D and N extrapolations, as well as the systematic error from choosing the fitting interval, while the second error is the difference between results from cubic and quartic fits in 1/ √ x. Note that the latter error . The dashed line corresponds to the fit (4.2), with the resulting value indicated on the axis. The lowest row corresponds to the ground state energy density, E0 2xN , as a function of 1/N , and is fitted according to (4.1). Notice that the lattice sizes are already close to the thermodynamic limit. The error bars are extracted from the extrapolation in D, shown in figure 4. is important only in the case of the ground state energy -for the scalar and vector mass gaps it is negligible with respect to the former error, as the data are very well described by polynomials only quadratic in 1/ √ x. Our result for the ground state energy in the massive cases is compatible with the massless one, and independent of m/g, as expected [35].
The results for the vector binding energies, M V g := ω 1 − 2m/g, are shown in table 1 for each value of m/g studied, together with the DMRG estimates [20] for comparison. In the massless case, the exact value is also displayed.

JHEP11(2013)158
Vector binding energy m/g MPS with OBC DMRG result [20] exact 0 0.56421(9) 0.5642 (2)   For the scalar, the most precise results found in the literature for the binding energy, M S g := ω 2 −2m/g in the massive case, correspond to the strong coupling expansion [36]. We show them together with our best fit and the exact value for the massless case in table 2.
One of the factors explaining the lower precision attained in the scalar calculation as compared to the vector results is the longer time required to reach the scalar state with the MPS algorithm for excited states. As discussed in section 3 the computational cost scales like the square of the required number of levels. Since in many instances the first scalar candidate appears at level 7 − 8 or above, this represents a cost over 50 times larger than in the case of the ground state. We have thus opted for a trade-off between precision and efficiency, and applied for the scalar a less demanding convergence criterion than for the vector or the ground state, what also translates in less precision in the final results. For the same reason, the largest bond dimensions are in some cases not available in the scalar calculations, which leads to somewhat less precise finite-size scaling.

Discussion
We have successfully employed open boundary MPS to compute the ground state and several excited levels of the lattice Schwinger model, using a staggered formulation, and for various values of the fermion mass. Although in the physical subspace in which we work the model contains long range interactions which do not decay with the distance, we have found that MPS produce a very good approximation not only to the ground state, but also to higher excited levels, and we are able to reach precisions comparable to those available from other techniques, or in some cases even slightly better. Additionally, the precision Excluding small values x < 50 produces a better fit, as shown by the solid line, for the interval x ∈ [100, 600]. Our final estimate (see appendix A) for the massless case is ω 0 = −0.318338 (22)(24), and for all other masses it is compatible with this value.
we reach is not extremely sensitive to the value of the mass, what further points to the usefulness of TN techniques for the non-perturbative parts of the parameter space.
In order to determine the mass gaps, we have also extended the MPS tools to approximate excited states, and we have proposed a modified algorithm that is efficient and allows us to approximate more than ten excited states in chains of hundreds of sites. These techniques provide an ansatz for the targeted states, so that not only the energies, but also other observables can be precisely determined [50].
Our results validate the applicability of tensor network techniques, in particular MPS, for lattice gauge theory problems. We have obtained very precise results, even though the open boundary MPS cannot represent states with well-defined momentum, a problem that a TI ansatz might likely overcome.  The long term goal would be to generalize these studies to higher dimensional systems, getting in this way closer to the target of lattice QCD studies. The MPS ansatz employed in this work is one-dimensional, and its applicability to higher dimensional problems is limited, a fact well understood in terms of its entanglement content. Nevertheless, the entanglement structure in the MPS has a natural generalization to two and higher dimensions in the PEPS ansatz [13,40]. Algorithms exist for the ground state search and for simulating the evolution also in the case of PEPS, and the ansatz can directly be applied to fermionic systems [51][52][53]. Although the higher computational costs limit the nowadays feasible system sizes and bond dimensions, PEPS are good candidates to study higher dimensional lattice gauge problems. Our MPS study is a first step to assess the feasibility of TNS to describe the relevant physics in this kind of problems.  One of the main advantages of the MPS and other TNS methods is that they easily allow us to attack other problems which are more complicated for standard lattice techniques. In particular, it is straightforward to include a chemical potential term in the Hamiltonian, which is interesting in many-flavour Schwinger models. Additionally, we could perform a similar study in a more complicated lattice theory, in particular with a non-Abelian symmetry. Thermal equilibrium properties can be easily studied with TNS, and it is also possible to simulate real time evolution, which opens the door to out-of-equilibrium questions. Thus even if tensor networks are presently not competitive with standard Monte Carlo techniques for 4D quantum field theories, they might allow us to address problems that are not amenable for customary lattice field theory techniques.

JHEP11(2013)158
Acknowledgments This work was partially funded by the EU through SIQS grant (FP7 600645). K.C. was supported by Foundation for Polish Science fellowship "Kolumb". This work was supported in part by the DFG Sonderforschungsbereich/Transregio SFB/TR9. K.J. was supported in part by the Cyprus Research Promotion Foundation under contract ΠPOΣEΛKYΣH/EMΠEIPOΣ/0311/16.

A Continuum limit and error analysis
As we have described in the main body of the text, we have a set of data corresponding to a given bond dimension D, lattice size N and coupling x. Having performed extrapolations to infinite D and N , we are left with the dependence of our observables on x and we want to take the continuum limit x → ∞. The leading order cut-off effects are O(1/ √ x), but we are sensitive also to higher-order corrections (O(1/x), O(1/x 3/2 ), . . .). Therefore, in general, our continuum extrapolated result depends on the chosen fitting interval in x. To take this systematic uncertainty into account, we perform the following error analysis procedure, similar to those employed in data analysis for Lattice QCD (see e.g. refs. [54,55]). Let us assume that, for a fixed value of m/g, we have a set of data for the chosen observable, O: for n values of x, such that O i is the observable evaluated at a finite lattice spacing corresponding to coupling x i , and ∆O i is its respective error (originating both from extrapolations to infinite D and N ). We choose some minimal number of points (n min ) and we fit every possible subset of p consecutive data points, with p ≥ n min , to a given polynomial function g(x). Each such fit gives us an estimate of the continuum value of the chosen observable, O c -we denote these continuum values by O c α , where α = 1, . . . , M fits . Each fit has its respective value of χ 2 , denoted by χ 2 α : where the sum runs over all data points entering the fit labelled by α. Our preferred value for the observable in the continuum limit is then defined as the median of the distribution of estimators O c α weighted by exp(−χ 2 α /N d.o.f. ), where N d.o.f. is the number of degrees of freedom.

(A.2)
Our preferred value for the considered observable in the continuum limit is then defined as the value of O α for which α is the smallest number that satisfies f α ≥ 0. f α ≥ 0.1585. This definition of the error is motivated by the fact that in the limit of an infinite number of fits, the weighted distribution will be Gaussian and this definition of the error will correspond to the width of such Gaussian distribution.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.