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.


I. INTRODUCTION
In the last years, methods based on tensor network states (TNS) have revealed themselves as very promising tools for the numerical study of strongly correlated quantum many-body systems. They are ansätze for the quantum state, characterized by their entanglement content and well suited for lattice systems. The paradigmatic family of TNS is that of matrix product states (MPS) [1][2][3][4][5][6], which underlies the well-known density matrix renormalization group algorithm (DMRG) [7,8]. The enormous success of DMRG for the study of one dimensional condensed matter problems has been better understood within the framework of TNS (see e.g. [9]). This has also enabled different extensions of the original algorithm, such as time-dependence [10][11][12], so that MPS/DMRG methods constitute nowadays a quasi exact method for the study of ground states, low lying excitations and thermal equilibrium properties of quantum spin chains far beyond the reach of exact diagonalization. Being free from the sign problem which plagues quantum Monte Carlo (QMC) methods, higher dimensional TNS [13][14][15] are seen as powerful candidates for the numerical exploration of long standing strongly correlated electron problems. This success has motivated the application of TNS techniques to further problems, and a relatively new field of application has been found in quantum field theories. Conformal Field Theory inspires a generalization of MPS [16] useful for critical models. Furthermore, specific generalizations of TNS exist [17][18][19] which are suitable for non-relativistic and relativistic QFT in the continuum. The discrete versions, on the other hand, are adequate for lattice field theories, and can in particular be applied to models relevant to high energy physics problems. This route was first explored using the original DMRG formulation [20,21]. The deeper understanding of the technique achieved thanks to TNS has increased the power of these methods and new results are now available for critical QFT [22,23].
A particularly interesting study case is that of the Schwinger model [24,25], or QED in one spatial dimension, the simplest gauge theory, which nevertheless exhibits features in common with more complex models (QCD) such as confinement, a non-trivial vacuum and asymptotic freedom, and has been adopted as a benchmark model where to explore lattice gauge theory techniques (see e.g. [26][27][28][29][30] 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 [31,32]. 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. 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 II we briefly introduce the Schwinger model, and review its formulation as a spin Hamiltonian. Section III presents the MPS formalism, with particular emphasis on the new techniques used to overcome the challenge posed by this particular kind of problems. In section IV we present our numerical results and conclude in section V with a discussion and outlook.

II. THE MODEL
The massive Schwinger model [24], 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 µ , g is 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 [24], 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 [25,33].
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 mass the best results were obtained with DMRG [20], while strong coupling expansion (SCE) got the most accurate prediction for the scalar mass [34,35] and the most precise values in the massless case [36].
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 [25].
The model can be formulated on a lattice. Here we focus on the Kogut-Susskind staggered formulation [37], 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 (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 (4) determines the electric field up to a constant α which can be added to L n and represents the background electric field.
Using a Jordan-Wigner transformation, φ n = k<n (iσ z k )σ − n , where σ ± = σ x ± iσ y , model (3) can be mapped to a spin Hamiltonian [31], Gauss' law reads then L n − L n−1 = 1 2 [σ z n + (−1) n ], and can be used to eliminate the gauge degrees of freedom, leaving the Hamiltonian [38] H =x 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.

III. METHOD
In this work we use the MPS ansatz to approximate the ground and lowest excited states of the Hamiltonian (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 [40], 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 [41] and recently extended in [42]. 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 [43].
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 [44,45]. However, finite-size effects are much larger for OBC and therefore simulation of larger chains may be needed to reach the thermodynamic limit reliably.
To determine the first excited state, after having found the ground state of the system, |Ψ 0 , we can construct the projector onto the orthogonal subspace, Π 0 = 1 − |Ψ 0 Ψ 0 |. The first excitation corresponds then to the state that minimizes the energy of the projected Hamiltonian, given that E 1 < 0, what we can always ensure by adding an appropriate constant to H.
This minimization corresponds to finding the ground state of the effective Hamiltonian 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 (9), which can be constructed from the set of all previously computed MPS levels {|Ψ k }. The expression to be minimized at each step of the MPS iteration is then 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) [46]. 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 computational cost 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 [50], 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 [47] 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 (Fig. 2).
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 [25], 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 respectively containing the vector, and the ground state with the scalar [48]. 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 [38] this was done 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 [51]. 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 2.

IV. RESULTS
We have applied the methods described in section III to the Hamiltonian (6) to determine the ground state, and the vector and scalar mass gaps for different 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.
The continuum limit corresponds to the lattice spacing a → 0, or x → ∞. The goal is then to compute the ground state energy density, and the vector and scalar mass gaps and estimate the error as the difference between this value and the one for the immediately smaller D). In figure 3 we show the results of this step for one particular case.
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 finite-size to extract the ground state energy density, E 0 /(2xN), and the mass gaps, ( 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 ) [38]. Hence the bulk limit is extracted by fitting the ground state energy density as and the energy gaps as We have run these fits for each pair of values (m/g, x), as illustrated in figure 4 for a particular case.
Finally, the values in the continuum limit can be extracted from the ones in the thermodynamic limit by fitting them, for each value of the mass, to a polynomial in 1 √ x . For these fits, we include only those values of x for which the thermodynamic limit can be extracted accurately (i.e. the corresponding level for the various system sizes is converged to a large enough D). The continuum limit extrapolations are shown in figures 5 (ground state), 6 (vector) and 7 (scalar). The interval [x min , x max ] and the function used for the fit affect the final result. We have thus used a quadratic function in 1 x and performed the fit over various intervals. We choose as our final estimation the value obtained from the best fit, and estimate the error from the difference with respect to another interval (see figure   captions).
Finally, we obtain for the ground state energy density in the massless case the value ω 0 = −0.31837 (7), to be compared to the exact value, −1/π ≈ −0.3183099. Our result for the ground state energy in the massive cases is independent of m/g, as expected [34]. The results for the vector binding energies, M V g := ω 1 − 2m/g, are shown in the following table for each value of m/g studied, with the DMRG estimates [20] shown for comparison. In the massless case, the exact value is also displayed.  [49].
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.
One of the main advantages of the MPS methods used in this work is that they easily allow us to attack other problems which are more complicated for standard lattice techniques. In particular, we can trivially include a chemical potential term in the Hamiltonian, and we can easily study thermal equilibrium properties. It is also possible to simulate real time evolution, which opens the door to out-of-equilibrium questions. Additionally, we could perform a similar study in a more complicated lattice theory, in particular with a non-Abelian symmetry, getting in this way closer to the target of lattice QCD studies. Thus, although 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.     for the scalar (blue) and vector (green) mass gaps, x . The dashed line corresponds to the fit (12), with the resulting value indicated on the axis. The lowest row corresponds to the ground state energy density, E 0 2xN , as a function of 1/N , and is fitted according to (11). 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 3.  shows the exact solution ω 2 (m/g = 0) = 2/ √ π ≈ 1.12838. Fits for the whole range of x (cubic) and for a smaller interval (x ∈ [5, 300]) are shown (respectively dashed and solid lines). As in figure 6, some of the large x estimations have large errors, because they require solving larger systems, for some of which the maximum bond dimension reached was not large enough. Again, we keep as estimate the fit for x ∈ [15,300], excluding the very small x, which do not follow the asymptotic behavior, and the larger values of x which have higher errors, and get a rough estimation of the error by comparing to [30,200] (except in the case m/g = 0.5, when due to the larger errors of relatively large x, we compare to the fit in [15,200]).