Topological vacuum structure of the Schwinger model with matrix product states

We numerically study the single-flavor Schwinger model with a topological θ -term, which is practically inaccessible by standard lattice Monte Carlo simulations due to the sign problem. By using numerical methods based on tensor networks, especially the one-dimensional matrix product states, we explore the nontrivial θ -dependence of several lattice and continuum quantities in the Hamiltonian formulation. In particular, we compute the ground-state energy, the electric field, the chiral fermion condensate, and the topological vacuum susceptibility for positive, zero, and even negative fermion mass. In the chiral limit, we demonstrate that the continuum model becomes independent of the vacuum angle θ , thus respecting CP invariance, while lattice artifacts still depend on θ . We also confirm that negative masses can be mapped to positive masses by shifting θ → θ þ π due to the axial anomaly in the continuum, while lattice artifacts nontrivially distort this mapping. This mass regime is particularly interesting for the ( 3 þ 1 )-dimensional QCD analog of the Schwinger model, the sign problem of which requires the development and testing of new numerical techniques beyond the conventional Monte Carlo approach. DOI:


I. INTRODUCTION
QED in 1 þ 1 dimensions, also known as the Schwinger model [1], is analytically solvable in the masslessfermion limit and exhibits many properties similar to QCD: confinement, chiral symmetry breaking, a Uð1Þ A quantum anomaly, and a topologically nontrivial vacuum leading to a θ-term. Therefore, the lattice-regularized version of the Schwinger model has been adopted as a benchmark model for developing and testing new numerical techniques. These comprise algorithms to tackle the sign problem, new ideas for Markov chain Monte Carlo (MCMC) investigations, and tensor network approaches (see, e.g., Refs. [2][3][4], respectively, and references therein).
Already since the seminal work by Coleman and collaborators [5,6], the role of topology in the Schwinger model has been extensively discussed. The nontrivial topological vacuum structure gives rise to a θ-dependent electric background field [5,6], which linearly depends on the fermion mass and gets completely screened in the zero-mass limit [7]. Similarly, the ground-state energy density, the chiral fermion condensate, and the topological vacuum susceptibility of the Schwinger model are nontrivially dependent on θ and the mass parameter. In particular, the topological vacuum susceptibility, which in the QCD analog measures the strength of CP violation of the theory, vanishes in the massless limit.
One well-studied regime in the parameter space of the Schwinger model is the second-order phase transition that occurs at θ ¼ π [6] and a fermion mass of m ≈ 0.33g [8,9], where g is the dimensionful coupling. A less studied regime is a vanishing or even negative fermion mass, which we consider in the current paper.
The motivation for looking at this parameter range is twofold. First, it is expected that a negative mass can be trivially mapped to a positive mass by shifting θ → θ þ π, which is given in the continuum due to the axial quantum anomaly [10][11][12]. However, as we demonstrate, lattice artifacts distort this mapping, so that negative-mass results only reproduce positive-mass results with θ → θ þ π after the continuum extrapolation. The negative-mass scenario is particularly interesting in the many-flavor case, where it can give rise to the CP-violating Dashen phase and pion condensation [13].
Second, the zero-mass regime is motivated by the possibility of having a vanishing up-quark mass in the QCD analog, which has been studied for several decades because it provides a potential solution of the long-standing strong CP problem [14][15][16][17]. The key point of this proposal is that topological effects can give rise to an effective up-quark mass that does not spoil the solution to the strong CP problem but nevertheless appears in the chiral Lagrangian. There is currently no evidence that this effective mass term is large enough to make the proposal phenomenologically viable (see Ref. [18] for a review), but, if feasible, it would be an elegant solution of the strong CP problem. Both this solution and the alternative axion solution of the strong CP problem [19][20][21] render the theory independent of the vacuum angle θ in the continuum. However, lattice artifacts still depend on θ, as we show for the single-flavor Schwinger model. Since all the above questions are of nonperturbative nature, it would be natural to address them through numerical lattice calculations, e.g., by the successful MCMC method. However, the Schwinger model has a sign problem when a topological θ-term is added to the action. This renders the MCMC approach inapplicable, at least when the value of the vacuum angle θ becomes too large.
In addition, as discussed above, we are also interested in going to zero and even negative fermion mass. In this case, the lattice Dirac operator D will develop zero modes. This is problematic for standard MCMC methods which use the operator D † D to have a real and positive action and then integrate out the fermions. The resulting determinant of D † D, is estimated stochastically (see, e.g., Ref. [22]) with appropriately chosen bosonic fields Φ † and Φ. Thus, if a zero mode appears in D † D, the integral in Eq. (1) is ill defined. If we also add a topological θ-term to study the physics questions described above, we encounter a double sign problem, thus ruling out a treatment with MCMC. A possible way out are tensor network techniques and, since we consider the (1 þ 1)-dimensional Schwinger model, in particular the matrix product state (MPS) approach. Investigations of gauge theories in 1 þ 1 dimensions, especially of the Schwinger model, have progressed substantially over the last years. There have been several works concentrating on the spectrum of the Schwinger model using MPS [4,9,[23][24][25][26][27][28][29][30][31][32]. The model has also been studied at nonzero temperature [4,[33][34][35][36][37], nonzero chemical potential [30,31,38] and for real-time problems [27,39]. In addition, quantum link models [40,41], Z n -QED models [42][43][44], and non-Abelian gauge models have been explored with the MPS approach [45][46][47][48][49][50]. Besides MPS, tensor network renormalization techniques [51][52][53] also have been very successfully employed to study properties of gauge theories in 1 þ 1 dimensions and recently even in a simple (2 þ 1)-dimensional gauge theory [54].
An early work on the Schwinger model with a topological term using density matrix renormalization group methods can be found in Ref. [8] and a more recent one using MPS in Ref. [9]. However, in these papers the main interest has been to explore the phase transition at a value of the vacuum angle θ ¼ π. Reference [9] nicely demonstrated the breaking of the CP symmetry happening at a first-order phase transition.
In the present work, we are interested in a different regime of the Schwinger model with a topological θ-term, namely at positive and negative fermion masses close to zero, i.e., far away from the phase transition. We aim to compute the θ-dependence of the ground-state energy density, the electric field, the chiral fermion condensate, and the topological vacuum susceptibility in this regime, in particular for the CP-conserving case of a vanishing mass. This extended regime allows us to examine the full range of validity of the mass-perturbation theory computations in the Schwinger model [7], including its limitations on the lattice. Our work also serves as a proof of concept that the MPS approach works well even for negative masses, which has not been explored before. The examination of this regime is particularly important for the long-term goal of applying such new numerical techniques to the abovementioned unresolved issues of QCD.

A. The Schwinger model
The massive Schwinger model [1] describes (1 þ 1)dimensional quantum electrodynamics coupled to a single massive Dirac fermion, with the Lagrangian density Here, ψ denotes the two-component fermionic field with bare mass m and A μ is the Uð1Þ gauge field with coupling constant g and field strength . The θ-term in Eq. (2) with θ ∈ ½0; 2π is a total derivative and therefore does not affect the classical equations of motion, but it does affect the quantum spectrum. In terms of the dimensionless parameter m=g of the model, both the massless case m=g ¼ 0 and the free case m=g → ∞ can be solved analytically. Therefore, the limits of very small or very large masses can be studied within perturbation theory, while the intermediate regime requires a nonperturbative treatment. The Hamiltonian density of the massive single-flavor Schwinger model in the temporal gauge, A 0 ¼ 0, reads plus an irrelevant constant. The electric field, F ¼ − _ A 1 , is fixed by the Gauß constraint, ∂ 1 F ¼ gψγ 0 ψ, up to an integration constant of gθ=2π, which corresponds to an electric background field [6]. The θ-parameter can be shifted between the electric field and the fermion mass term by performing an anomalous axial rotation of the fermionic field (see Appendix A for details). The axial quantum anomaly is also the reason why the Hamiltonian in Eq. (3) contains a θ-term at all, even though this term can be stripped away on the classical level when the Hamiltonian is formulated in terms of the electric field [55,56].
The dependence of several observables on the constant electric background field gθ=2π was computed in the continuum Schwinger model for the limit m=g ≪ 1 with mass-perturbation theory [7]. The ground-state energy density in units of the coupling was found to be [57] E 0 ðm; θÞ where Σ ¼ ge γ =ð2π 3=2 Þ, γ is the Euler-Mascheroni constant, and μ 2 0 E þ ¼ −8.9139, μ 2 0 E − ¼ 9.7384 are numerical constants. Note that the topological cosine structure of the ground-state energy density appears analogously in QCD and axion physics, and plays an important role for axion phenomenology (see, e.g., Ref. [59]).
From the energy density (4), one can obtain the electric field density, which is (up to a factor) its derivative with respect to θ. In units of the coupling, it is given by Thus, in the presence of massive fermions, the constant electric background field density F ðθÞ ¼ gθ=2π gets partially screened due to vacuum polarization. In the massless case, the field is completely screened, which can be equivalently described by the elimination of the vacuum θ-angle by an axial fermion rotation. The second derivative of the energy density with respect to θ corresponds (up to a sign) to the topological vacuum susceptibility. In units of the coupling, it reads This quantity vanishes in the chiral limit where physics becomes θ independent, similar to the QCD case, where the topological vacuum susceptibility is a measure for CP violation. Note that the susceptibility in Eq. (6) as well as the ground-state energy in Eq. (4) and the electric field in Eq. (5) are invariant under the simultaneous shifts of θ → θ þ π and m → −m. This is because the θ-parameter can be rotated from the electric field into the fermion mass term, as explained above. Thus, a shift by Δθ ¼ π gets compensated for by a change of the mass sign, m expðiΔθÞ ¼ −m.
Finally, one can also compute the chiral fermion condensate C ¼ hψψi, which is given by the derivative of the energy density with respect to the bare fermion mass, and which is independent of the fermion mass in first order. The condensate transforms similarly to the fermion mass under an axial transformation, therefore the abovementioned shift of θ → θ þ π induces m → −m and Cðm; θÞ → −Cð−m; θ þ πÞ, as can be seen in Eq. (7). In the massless limit, the condensate still seems θ dependent at first sight, but this angular parameter becomes unphysical as it can be rotated away by said axial rotation. Equivalently, for m ¼ 0 the phase of the condensate can be absorbed by a shift in the Schwinger boson field and becomes unobservable. Also note that the condensate itself it not a physical quantity and only enters observable quantities when multiplied by the fermion mass. Thus, the model's θ-dependence still vanishes for m ¼ 0.

B. Lattice formulation
Our goal is to numerically compute the θ-dependence of the quantities in Eqs. (4)-(7) using the MPS approach, for positive, zero, and negative fermion masses. For our numerical simulations with MPS, we use a lattice formulation of the Schwinger Hamiltonian. To distinguish lattice from continuum quantities in our equations, we denote lattice quantities with roman letters as opposed to the calligraphic letters for continuum quantities.
A possible discretization of the continuum Schwinger Hamiltonian in Eq. (3) on a lattice with spacing a is given by the Kogut-Susskind Hamiltonian [60] In the expression above, ϕ n is a single-component fermionic field describing a fermion on site n, m is the bare fermion mass, and g is the coupling constant. The operators F n and ϑ n act on the gauge links in between the fermions, and F n gives the quantized electric flux on link n. They fulfill the commutation relation ½ϑ n ; F k ¼ iδ n;k and hence e iϑ n acts as rising operator for the electric flux. The angle ϑ n is restricted to ½0; 2π since we use a compact formulation. The Gauß constraint on the physical states translates to on the lattice, where Q n ¼ ϕ † n ϕ n − ð1 − ð−1Þ n Þ=2 is the staggered fermionic charge.
We choose to work with open boundary conditions, for which we can use Eq. (9) to integrate out the gauge fields [23,61]. After a residual gauge transformation, the dimensionless lattice Hamiltonian reads where we have defined the dimensionless constants x ≡ 1=ðagÞ 2 and μ ≡ ffiffi ffi x p m=g. As in the continuum case, the dimensionless integration constant θ=2π corresponds to a constant electric background field. Hence, we see that the model only has three independent parameters: the lattice spacing and the bare fermion mass, both in units of the coupling, and the vacuum angle θ.
We would like to make contact to the continuum prediction with the results extracted from the dimensionless lattice Hamiltonian. To this end, let E 0 ðm; θÞ be the groundstate energy of the dimensionless Hamiltonian W from Eq. (10), which is related to the dimensionful energy density E 0 ðm; θÞ from Eq. (4) by Thus, starting from E 0 ðm; θÞ, we find where we have used that the volume L in units of the coupling is given by Consequently, the lattice quantity that should correspond to the continuum energy density (4) is E 0 =2N. Note, however, that E 0 =2N is UV divergent [58]. In order to obtain a UV-finite quantity, we can simply subtract the result for a fixed value of θ ¼ θ 0 and look at ΔE 0 ðm; θÞ=2N ¼ ½E 0 ðm; θÞ − E 0 ðm; θ 0 Þ=2N. For small enough bare fermion mass, we expect that toward the continuum limit this UV-finite lattice ground-state energy density becomes approximately equal to the perturbative continuum prediction ΔE 0 ðm; θÞ [62], After extrapolating our numerical lattice data for the UV-finite ground-state energy density ΔE 0 ðm; θÞ to the continuum, we can numerically compute the derivatives to obtain the continuum electric field density F ðm; θÞ and the continuum topological vacuum susceptibility χðm; θÞ. The electric field density and the topological vacuum susceptibility are already UV finite without subtracting their values at θ 0 , which is why we denote them as F and χ instead of ΔF and Δχ, respectively. Alternatively, we can also directly measure the electric field per unit volume with MPS and numerically compute its derivative to get the continuum topological vacuum susceptibility. As ðF n þ θ=2πÞg approaches the electric field F ðxÞ in the continuum limit, we expect F n þ θ=2π to follow Eq. (5) for small fermion masses toward the continuum. Notice that we are using a staggered formulation; thus, in order to compensate for staggering effects, we average F n over the system and look at the quantity In addition, we also have direct access to the chiral condensate, which in our staggered formulation translates For nonvanishing bare fermion mass, this quantity is UV divergent [24,27,37]. Thus, we again subtract the value for θ 0 , ΔCðm; θÞ ¼ Cðm; θÞ − Cðm; θ 0 Þ, and expect to find for small fermion masses toward the continuum limit

C. Matrix product states
In order to obtain the ground state of the Hamiltonian in Eq. (10), we use the MPS ansatz. For a system with N sites and open boundary conditions, it reads The size D of the matrices, called the bond dimension of the MPS, determines the number of variational parameters in the ansatz and limits the amount of entanglement that can be present in the state (see Refs. [63][64][65] for detailed reviews).
Given a Hamiltonian, the MPS approximation for the ground state can be found in a standard manner by iteratively updating the tensors A i k k one after another while keeping the others fixed [66]. In each step, the optimal tensor is determined by finding the ground state of an effective Hamiltonian describing the interaction of the site with its environment. The ground-state wave function is obtained by repeating the updating procedure starting from the left boundary and sweeping back and forth until the relative change of the energy is below a certain tolerance η. After obtaining the MPS for the ground state, we can measure all kinds of (local) quantities such as the electric field and the chiral condensate.
For convenience in the simulations, we choose to translate the fermionic degrees of freedom in Eq. (10) to spins using a Jordan-Wigner transformation [23]. Although tensor networks and in particular MPS can deal with fermionic degrees of freedom with essentially no additional cost in the algorithm [67][68][69], this allows us to avoid dealing with anticommuting fermionic operators. In spin language, Eq. (10) reads where Q n ¼ ðσ z n þ ð−1Þ n Þ=2 is the staggered charge and σ AE ¼ ðσ x AE iσ y Þ and σ z are the usual Pauli matrices.
Although the Hamiltonian in Eq. (17) is nonlocal, we expect MPS to be a suitable ansatz to describe its lowenergy spectrum. For one, the original Hamiltonian (8) is local, and for bare fermion masses smaller than ðm=gÞ c ≈ 0.33 the model is gapped with a unique ground state. Moreover, as shown in Ref. [9], its low-energy states are characterized by small values of the electric field. Hence, the gauge links can be effectively considered as finite dimensional and the ground state can be efficiently described by a MPS [70]. Integrating the gauge degrees of freedom out can be seen as projecting the ground state of the original Hamiltonian locally on the gauge invariant subspace. Such a projection increases the bond dimension at maximum by a factor on the order of the effective dimension of the gauge links. As a result, the ground state of Hamiltonian (17) is expected to be well described by a MPS with small bond dimension, too.

III. RESULTS
We examine the θ-dependence of the ground-state energy density, the electric field, the chiral condensate, and the topological vacuum susceptibility for a wide range of fermion masses, m=g ∈ ½−0.07; 0.21, and lattice spacings corresponding to x ¼ 1=ðagÞ 2 ∈ ½80; 160. In order to systematically probe for finite-volume effects and to be able to extrapolate our results to the thermodynamic limit, we explore for each combination of ðθ; m=g; xÞ a large range of system sizes corresponding to volumes N= ffiffi ffi x p ∈ ½4.5; 45. In addition, we have another truncation effect due to the finite bond dimension present in our numerical simulations. This error can be controlled by repeating the calculation for every set of ðθ; m=g; x; NÞ for a range of D ∈ ½20; 140 and extrapolating to the limit D → ∞ (details about the extrapolation procedure can be found in Appendix B). In all our simulations, we stop as soon as the relative change in the ground-state energy is below η ¼ 10 −10 . Moreover, we focus on half a period of θ ∈ ½0; π, since all observable quantities are predicted to be (point) symmetric around θ ¼ π [see Eqs. (4)- (7) and exemplary data of a full period of θ ∈ ½0; 2π in Appendix C].

A. Ground-state energy
Our results for the UV-finite ground-state energy density are shown in Fig. 1, after subtracting the value for θ 0 ¼ 0 and extrapolating to the limit N → ∞. The figure contains our data for the largest and smallest lattice spacing as well as the continuum extrapolation. Note that the result for θ 0 ¼ 0, which is subtracted in ΔE 0 ðm; θÞ, is smaller than FIG. 1. UV-finite ground-state energy density as a function of the angle θ for (a) m=g The red dots represent the result obtained after extrapolating our finite-lattice data to the continuum which can be compared to the perturbative prediction from Eq. (13) (blue solid line). In all cases the error bars are smaller than the markers.
the results for θ > θ 0 . This is why the UV-finite groundstate energy density in Fig. 1 is positive for m > 0, while the UV-infinite expression in Eq. (4) is negative.
In general, we observe that we can reliably extrapolate to the thermodynamic limit and control our errors due to the finite bond dimension and system size (see Appendix B for details). Independent of the fermion mass, we see that finite-lattice effects are more pronounced around the extremal values of the energy density at θ ¼ π. For small masses, our data exhibit larger relative changes when going to smaller lattice spacings, whereas there is hardly any shift for the two largest masses m=g ¼ 0.14 and 0.21, especially for small θ.
With our finite-lattice data, we can extrapolate to the limit ag → 0 and estimate the continuum values. As Fig. 1 reveals, the error bars resulting from this extrapolation are negligible and we can obtain precise estimates for the continuum limit. Comparing our results to the prediction from mass-perturbation theory in Eq. (13), we observe excellent agreement for For the largest two masses, perturbation theory eventually breaks down and fails to describe our data, as expected. Even though m=g ¼ 0.21 is still much smaller than unity, the perturbative results become less reliable than our numerical computations in this regime. This is illustrated by comparing the perturbative prediction of the critical mass of the phase transition, ðm=gÞ c ≈ 0.18 (based on Ref. [71]), with the nonperturbative value ðm=gÞ c ≈ 0.18 ≪ 0.33. Although the perturbative calculations in Ref. [7] are still qualitatively correct and approximately follow our numerical data even in the large-mass parameter regime [see Figs. 1(d) and 1(e)], we can only obtain quantitatively precise results with numerical techniques beyond perturbation theory.
In particular, it is interesting to look at our continuum data for vanishing bare fermion mass. As predicted by the perturbative result in Eq. (13), we indeed observe that the energy density becomes independent of θ once the extrapolation to the continuum is performed. This elimination of the θ-parameter in the chiral limit is only given in the continuum and does not apply to finite-lattice spacings, as one can see in Fig. 1(b).
Moreover, it is instructive to compare our results for m=g ¼ −0.07 and 0.07. In the continuum, a negative bare fermion mass can be mapped to the same positive mass value by shifting the θ-angle by π [see Eq. (13)]. Hence, a negative mass yields the same ground-state energy for θ ∈ ½0; π as the corresponding positive mass for θ ∈ ½π; 2π. This can be seen when comparing Figs. 1(a) and 1(c), keeping in mind that the mapping requires not only the shift θ → θ þ π but also the shift θ 0 → θ 0 þ π, i.e., the UV-finite quantity is obtained by subtracting the value at θ 0 ¼ π instead of θ 0 ¼ 0.
Equivalently, the continuum data for m < 0 can be mapped to the continuum data for m > 0 by reflecting the former with respect to the x axis. Even though this mapping works well in the continuum, Fig. 2 reveals that it is distorted for our finite-lattice data. This is expected due to the difficulty to establish chiral symmetry and the corresponding axial anomaly on the lattice [72].
The artifacts from the finite-lattice spacing enter with opposite sign, breaking the reflection symmetry for negative and positive masses. Moreover, for a fixed value of the lattice spacing, the deviations from the continuum result not only differ in sign, but also their absolute value , we see that these differences disappear when extrapolating to the limit of vanishing lattice spacing. Thus, the reflection symmetry is restored and our continuum data is in excellent agreement with the perturbative prediction for small masses.

B. Electric field
With our MPS approach, we can directly measure the electric field in the ground state. Performing the same extrapolation procedure as for the energy density, we obtain the data shown in Fig. 3. Again, the errors resulting from the extrapolation in system size and bond dimension are negligible, and finite-lattice effects are more pronounced for smaller values of m=g. The lattice effects become stronger for θ ¼ π=4 due to the sine dependence of the electric field [see Eq. (5)], in contrast to the ground-state energy in Fig. 1 whose cosine-dependent lattice effects become more distinct for θ ¼ π [see Eq. (4)].
Just as before, we can estimate the continuum value of the electric field by extrapolating our finite-lattice results to the limit of vanishing lattice spacing. In general, we find that we can reliably estimate the continuum limit, while the errors are slightly larger compared to the energy density. In particular, for the largest two masses we observe enhanced error bars for a θ-angle of π, which is likely caused by the fact that the continuum extrapolation becomes increasingly challenging as we approach the phase transition at ðm=gÞ c ≈ 0.33 and θ ¼ π.
When comparing our continuum results for the electric field to the predictions from mass-perturbation theory [see Eq. (5)], we observe a similar picture as for the ground-state energy density. For bare fermion masses m=g ≤ 0.07 [see Figs. 3(a)-3(c)], the perturbative prediction is in excellent agreement with our continuum results. In particular, for m=g ¼ 0 our data are compatible with zero, independent of θ. Hence, our results confirm that for vanishing mass, the background field gets screened completely and the total electric field vanishes. Moreover, our data for m=g ¼ −0.07 and m=g ¼ 0.07 can be mapped in the continuum by reflection with respect to the x axis, while finite-lattice artifacts distort this mapping with opposite sign. Moving on to larger fermion masses, perturbation theory eventually breaks down and Eq. (5) does not reproduce the behavior of our data for m=g ≥ 0.14. In contrast to the energy density, the perturbative result does not even qualitatively describe our numerical data for the largest mass, m=g ¼ 0.21. As pointed out before, this mass regime cannot be accurately captured by perturbation theory, which manifests itself in predicting a wrong critical mass of the phase transition.
As a cross-check, we can also obtain data for the electric field density by numerically computing the derivative of our results for the continuum energy density, which we show in Fig. 3 for comparison. In general, the values for the electric field obtained in this way are in good agreement with those from extrapolating the direct measurement of the electric field to the continuum. Although numerically computing the derivative enhances the errors by a factor proportional to 1=Δθ, the results from our data for the energy density are typically more precise than the extrapolated electric field values. The reason for this is that the electric field is in general more sensitive to finite-volume and finite-lattice effects (see Appendix B), resulting in larger error bars in the extrapolated values. For large fermion masses, the errors increase in both cases around θ ≈ π, thus indicating that we get closer to the critical value ðm=gÞ c .

C. Chiral condensate
The MPS approach also gives us access to the chiral condensate in the ground state. Performing the same extrapolation procedure as for the ground-state energy density and the electric field, we obtain the results in Fig. 4, where we have again subtracted the value for θ 0 ¼ 0. As before, the result for θ 0 ¼ 0, which is subtracted in ΔC 0 ðm; θÞ, is smaller than the results for θ > θ 0 . This is why the UV-finite chiral condensate in Fig. 4 is positive for m > 0, while the UV-infinite expression in Eq. (7) is negative.
Compared to the ground-state energy density and the electric field, the chiral condensate is less susceptible to finite-lattice effects and there is hardly any difference between results for our coarsest and finest lattice spacing. As a result, the error bars from extrapolating our finitelattice data to the limit ag → 0 are essentially negligible.
Notice that the chiral condensate corresponds to the derivative of the energy density with respect to the bare fermion mass [see Eq. (7)]. Thus, in leading order, massperturbation theory predicts a behavior independent of the parameter m=g. Looking at our results for small masses in Figs. 4(a)-4(c), we see that this indeed is the case, and Eq. (15) is in excellent agreement with our data. At next order in perturbation theory, we expect finitemass effects to become relevant, which becomes particularly interesting when comparing our data for m=g ¼ −0.07 and 0.07. Since the chiral condensate transforms similarly to the fermion mass under an axial rotation, a shift of θ → θ þ π does not only induce m → −m but also Cðm; θÞ → −Cð−m; θ þ πÞ [see Eq. (15)]. Hence, changing the sign of the condensate for m < 0 in the range θ ∈ ½0; π reproduces the corresponding positive condensate for m > 0 in the range θ ∈ ½π; 2π. This can be seen in Figs. 4(a) and 4(c), keeping in mind that the shifted UVfinite quantity is obtained by subtracting the value at θ 0 ¼ π instead of θ 0 ¼ 0.
In the massless limit, our data for the chiral condensate are still θ dependent [see Fig. 4(b)], but θ becomes an unphysical parameter as it can be rotated away by the above-mentioned axial rotation (see Appendix A). For larger values of m=g, the data only change moderately [see Figs. 4(d) and 4(e)], whereas perturbation theory erroneously predicts new qualitative features that do not occur in our numerical data, such as a dip around θ=2π ≈ 0.19. Thus, similar to the electric field, the prediction by mass-perturbation theory for the chiral condensate in Eq. (7) breaks down for large masses m=g ≥ 0.14 as we approach the phase transition.

D. Topological vacuum susceptibility
Although we cannot directly measure the topological vacuum susceptibility in our Hamiltonian framework, we can obtain this quantity by either numerically differentiating our continuum estimates for the electric field or by computing the second derivative from our data for the ground-state energy [see Eq. (6) and Appendix B for details]. Figure 5 shows our results for both approaches.
In general, both methods give remarkably consistent results. The data obtained from the second derivative of the energy density have noticeably smaller error bars, except for m=g ≥ 0.14 and θ ≈ π. Already the electric field and the ground-state energy density showed that the perturbative result breaks down for our largest two values of m=g, and For small masses, our data are again in excellent agreement with the predictions from mass-perturbation theory, as Figs. 5(a)-5(c) reveal. In particular, Fig. 5(b) shows that for vanishing fermion mass our data are compatible with χ top =g ¼ 0, once more indicating that the background field gets completely screened in that case and θ is not a physical parameter. This becomes even more apparent in Fig. 6, where we show our data for the susceptibility as a function of the bare fermion mass for various values of θ [74]. The figure clearly shows that for vanishing bare fermion mass, the different curves intersect at χ top =g ≈ 0, which demonstrates the CP invariance of the Schwinger model for m=g ¼ 0. Just as in QCD, where the topological vacuum susceptibility is a measure of CP violation, the presence of a massless fermion allows the θ-angle to be rotated away by an axial fermion rotation, thus CP is preserved. The same rotation also maps our results for negative and positive masses when shifting θ → θ þ π [see Figs. 5(a) and 5(c)], as we already observed for the ground-state energy and the electric field. Note that Fig. 6 does not reveal this mapping, since the susceptibility is only shown for small values of θ ≪ π.
Finally, we point out that the topological susceptibility at θ ¼ 0 (corresponding to the blue markers in Fig. 6) becomes negative for negative fermion masses and is expected to diverge to negative infinity at m=g ≈ −0.33 (see Ref. [75] for a discussion of a similar effect in two-flavor QCD). This is because the well-known phase transition at m=g ≈ 0.33 and θ ¼ π is equivalent to a phase transition at m=g ≈ −0.33 and θ ¼ 0, due to the abovementioned mapping. Thus, the diverging susceptibility at θ ¼ 0 is associated with the diverging correlation length as one approaches the critical point. This phase transition would be nontrivial to study in the two-flavor Schwinger model with two masses of opposite sign, giving rise to the CP-violating Dashen phase [13].

IV. CONCLUSION
In this paper, we systematically explore the topological vacuum structure of the Schwinger model with a θ-term. In particular, we study the θ-dependence of the ground-state energy density, the electric field, the chiral condensate, and the topological vacuum susceptibility at zero and negative fermion mass. This mass regime is especially interesting for models with a sign problem, which require the development and testing of new numerical techniques beyond the conventional MCMC approach. The prime example would be (3 þ 1)-dimensional QCD, where the zero-mass case was proposed as a possible solution to the strong CP problem.
Addressing the Hamiltonian lattice formulation of the Schwinger model with numerical methods based on MPS, we show that we can reliably compute the ground state of the model in a controlled manner with small errors. While at very small masses we find excellent agreement with massperturbation theory, thus scrutinizing our approach, we also demonstrate the limitations of perturbative methods.
Our results provide us with a comprehensive picture of the topological vacuum structure of QED in 1 þ 1 dimensions. For small masses, the θ-dependence follows the perturbative analytical calculation from Ref. [7]. In the chiral limit, our data confirm that the model becomes CP invariant as the topological vacuum susceptibility vanishes and the θ-dependent electric background field gets screened due to vacuum polarization. Thus, for m=g ¼ 0, the θ-parameter that labels the topologically nontrivial vacua of the Schwinger model becomes an unphysical parameter, just as in the (3 þ 1)-dimensional analog of QCD with a massless up quark [14][15][16][17] or with an axion [19][20][21]. As we go to larger masses, the perturbative prediction eventually breaks down and especially the chiral condensate deviates significantly from it.
Comparing our finite-lattice and continuum data, we find that lattice artifacts reintroduce the θ-dependence of the observables in the massless limit. Moreover, in the massive regime, lattice artifacts enter inversely and with different strengths for opposite mass sign, which renders the negative-mass regime nontrivial on the lattice. In the continuum, negative masses can be trivially mapped to positive masses by shifting θ → θ þ π due to the quantum anomaly [10][11][12], which gets confirmed by our data. Our results demonstrate that MPS work well even for negative fermion masses, which has not been explored before and becomes particularly relevant in the many-flavor case, where a negative mass can generate a second-order phase transition to the CP-violating Dashen phase [13].
Regarding the topological vacuum angle θ, there are several interesting aspects that could be studied in the future. Generalizing our MPS setup for the Schwinger model to multiple flavors would be straightforward [30,31], allowing us to study the above-mentioned Dashen phase [13]. Moreover, our MPS approach is not limited to the Abelian case [4,45,[47][48][49][50] and could offer the possibility to explore non-Abelian gauge models in the presence of a topological θ-term. In this Appendix, we explain the origin of the θ-term in the Hamiltonian density of the continuum Schwinger model (3). On the classical level, the θ-term can be stripped away when the Hamiltonian is formulated in terms of the electric field [55,56]. However, Coleman argued on physical grounds that the θ-term can be understood as a background electric field in the quantized theory [6]. A rigorous derivation was done in Ref. [76], which also showed there are several technical subtleties. To make the paper selfcontained, we here present this computation, which gets simplified in the bosonized formalism and addresses all technical challenges arising in the Hamiltonian formulation. Extensions to higher dimensions would be possible following Ref. [77].

ACKNOWLEDGMENTS
The bosozined version of the continuum Schwinger Hamiltonian in Eq. (3) reads plus an irrelevant constant. Here, m ¼ g= ffiffiffi π p is the Schwinger mass and ΦðpÞ and ΠðpÞ are the bosonic operators that satisfy canonical commutation relations and hermeticity properties. In the following subsections, we will discuss the properties of these bosonic operators for p ≠ 0 and p ¼ 0 and compute the (non)conservation laws of the corresponding bosonic currents. Finally, we will use the anomaly equation (A18) obtained with the bosonized Hamiltonian (A1) to demonstrate that a constant shift in the electric field of the fermionic Schwinger Hamiltonian is equivalent to an axial rotation of the massive fermionic field.

Bosonic operators for nonzero momentum
The bosonic operators ΦðpÞ and ΠðpÞ of the Hamiltonian density (A1) are defined for p ≠ 0 as in terms of the bosonic chiral-charge density operators Here, the index α ¼ 1, 2 denotes the left-handed (α ¼ 1) and right-handed (α ¼ 2) states, p; k ∈ Z are the integer momenta, and the momentum-space Fermi operators a α;k are related to the real-space Fermi operators ψ α ðxÞ via ψ 1 ðxÞ The bosonic currents for p ≠ 0 are given by In the Coulomb gauge, ∂ x A x ¼ 0, the bosonic operator ΦðpÞ (A2) can be expressed in terms of the longitudinal part of the electric field, by using Gauß's law, and applying the relations ∂ x ¼ −ip and ∂ x A t ¼ −F long x ðpÞ to the left-hand side of Eq. (A9) and the current-operator relation (A6) to the right-hand side. Since the total electric charge is zero, the longitudinal part of the electric field F long x ðpÞ has no p ¼ 0 component, which is the quantization constraint in the Coulomb gauge.

Bosonic operators for zero momentum
The zero-momentum scalar field operators Φð0Þ and Πð0Þ read where VðA x Þ is the potential and F tr x ð0Þ is the transverse part of the electric field which is classically given by ∂ t A x . The transverse part of the electric field (A12) only contains the constant p ¼ 0 component, i.e., F tr x ð0Þ ¼ 2πF tr x .

(Non)conservation laws of bosonic currents
The Hamiltonian density (A1) yields the following equations of motion for the operators ΦðpÞ and ΠðpÞ: where we used d=dx ¼ −ip. This is the well-known conservation law for the vector current, where we used m ¼ g ffiffiffi π p and F x ¼ F long x ðpÞ þ F tr x ð0Þ. This is the well-known quantum anomaly equation for the axial current, We note that p 2 Φð0Þ ¼ −d 2 Φð0Þ=dx 2 vanishes on the lefthand side of Eq. (A17), therefore the zero-momentum contribution of ΦðpÞ only shows up in the nonderivative term on the right-hand side. For the nonderivative term, we have to add up both the p ≠ 0 and the p ¼ 0 contributions given by Eqs. (A8) and (A10).

Axial rotation of fermionic fields
Using the quantum anomaly equation (A18) obtained with the bosonized Hamiltonian, one can now consider the fermionic Hamiltonian (3) to demonstrate that an axial rotation of the massive fermionic field ψ, induces a constant shift in the electric field F x . The rotation (A19) by an infinitesimal angle θ ≪ 1 shifts the fermion mass term by and analogously the chiral fermion condensate by Since the fermion mass perturbatively corrects the divergence (A18) of the axial current j μ 5 ¼ψγ 5 γ μ ψ by the axial rotation of the fermionic field (A19) induces the following shift in the Hamiltonian density: where we used γ 5 γ 0 ¼ −γ 0 γ 5 . Thus, the angular parameter θ in the fermion mass term can be absorbed by a shift in the electric field, F x → F x þ θ=2π, and vice versa. In particular, for m ¼ 0 the θ-parameter becomes unphysical, since it can be rotated away without affecting any mass term, which corresponds to absorbing the phase of the chiral condensate hψψi in the Schwinger boson field.

APPENDIX B: DETAILS OF THE EXTRAPOLATION PROCEDURE
Here we give some details how we control the errors in our numerical simulations and how we extrapolate our data to obtain first the thermodynamic limit and finally the continuum limit. Since the extrapolation procedure has to be done independently for each combination of ðm=g; θÞ, we suppress all arguments in the following and just refer to E 0 =2N, F av , and C, meaning that we look at these quantities at a specific value of ðm=g; θÞ.
In a first step, for every combination of ðN; x; m=g; θÞ, we estimate the numerical error due to the finite matrix size in our MPS ansatz. To this end, we plot the observables O that we measure as a function of 1=D, and we use the three data points with the largest values for D to linearly extrapolate to the limit D → ∞ (see Fig. 7 for an example). We proceed in a standard manner: the central value is taken to be the mean value of our data point with the largest bond dimension, O D max , and the extrapolated value O D¼∞ . The error is estimated as half of the absolute value of their difference, δO D ¼ jO D max − O D¼∞ j=2. In general, we find that our bond dimensions are large enough to avoid noticeable truncation effects due to the finite matrix size. In addition, we have another error due to the finite convergence tolerance η in our simulations which results in δE 0;η ¼ ηE 0;D max for the ground-state energy and δO η ¼ ffiffi ffi η p O D max for other local observables [78]. The total error is then estimated as the square root of the sum of squares, After estimating the numerical errors due to the finite matrix size, we extrapolate our data for each combination of ðx; m=g; θÞ to the infinite-volume limit N → ∞, where we propagate our errors from the extrapolation in D. In general, we observe strong finite-volume effects for volumes N= ffiffi ffi x p < 15, thus for the extrapolation we only consider volumes larger than that. To estimate the infinite-volume limit, we fit our data to polynomials in 1=N up to degree 3 (see Fig. 8 for an example). For each polynomial, we try every fitting interval of consecutive data points, which contains at least two more data points than the degree of the polynomial. To obtain the central value, we choose the fit with the lowest value of χ 2 d:o:f: . In case we have several fits with χ 2 d:o:f: < 1, we choose the polynomial of smallest degree in 1=N that achieves this value. In most cases, we find that a linear fit in 1=N is enough to describe our data well. In addition to the error of the fitting coefficient, we estimate our systematic error. To this end, we compare our central value to the value obtained from the next best fit using the same degree polynomial or to the one obtained with the next highest order. The total error is then estimated, analogously to the extrapolation in D, as the square root of the sum of squares.
In a final step, we extrapolate our data to the continuum corresponding to ag → 0. To this end, we fit our finitelattice data again to polynomials in ag. In general, we observe that a linear function or sometimes even a constant is enough to describe our data well (see Fig. 9 for an example). We again propagate the errors from the extrapolation in N to estimate the final error of our data.
To obtain the electric field from our continuum estimates for the energy density, we follow Eq. (5) and compute the derivative numerically using second-order finite differences [79] For all the data we show, the distance between two different points is Δθ ¼ 0.025. The error of the electric field is then In all panels the red triangles correspond to the data used to extrapolate to the thermodynamic limit, the blue solid lines to the best linear fit, the orange dashed lines to the best quadratic fit, and the green dotted lines to the best cubic fit in 1=N.
estimated by propagating the error in the UV-finite energy density as a systematic error. The topological vacuum susceptibility can be calculated in a similar fashion, we can either numerically differentiate our results for the UV-finite energy density twice with respect to θ, or compute the derivative of our results for the electric field, where again the distance between two different angles is Δθ ¼ 0.025 and we propagate the errors in the UV-finite energy density and in the electric field values as systematic errors to obtain an error estimate for the topological susceptibility.

APPENDIX C: DATA FOR A FULL PERIOD OF θ
In the main text, we focused on the regime 0 ≤ θ ≤ π, since all quantities studied are (point) symmetric around π. Figure 10 shows an explicit example of a full period for m=g ¼ 0.21 after extrapolating to the thermodynamic limit. We restrict ourselves to a single lattice spacing corresponding to x ¼ 80, since this value is already very close to the continuum limit for such large masses, as demonstrated by the data in the main text. Looking at Fig. 10, we see that the data for the ground-state energy density, the electric field, and the chiral condensate are indeed symmetric around θ ¼ π, and we can obtain precise estimates throughout the entire period of θ ∈ ½0; 2π.
As theoretically predicted in Ref. [6] and numerically demonstrated in Refs. [8,9], the continuum model exhibits a first-order transition at θ ¼ π for bare fermion masses larger than the critical value ðm=gÞ c ≈ 0.33. This is accompanied by a spontaneous breaking of the CP symmetry, and the critical line ends in a second-order quantum phase transition exactly at ðm=gÞ c . Since our value for the largest bare fermion mass is still smaller than the critical one, we do not expect a transition to happen at θ ¼ π.
At first sight, our data for the chiral condensate in Fig. 10(c) give the impression that there could nevertheless be a transition, as we observe a sharp peak at that value of θ. However, due to the large slope of this quantity, our resolution in θ is limited. Taking a closer look at Fig. 10(b), we see that the electric field at θ ¼ π vanishes, thus indicating that the CP symmetry is not broken [6] and there is no phase transition. Our data for the energy density [ Fig. 10(a)] corroborates this picture. Although for θ ¼ π there is a noticeable peak, we do not observe a cusp in the data, thus giving another indication that for m=g ¼ 0.21 there is no transition. Nevertheless, the features in the energy density, the electric field, and the chiral condensate hint toward an upcoming phase transition at ðm=gÞ c ≈ 0.33.