Massive Schwinger model at finite $\theta$

Using the approach developed in [V. Azcoiti, G. Di Carlo, A. Galante, V. Laliena, \textit{Phys. Lett.} \textbf{B563}, (2003) 117], we are able to reconstruct the behavior of the massive 1-flavor Schwinger model with a $\theta$ term and a quantized topological charge. We calculate the full dependence of the order parameter with $\theta$. Our results at $\theta = \pi$ are compatible with Coleman's conjecture on the phase diagram of this model.


I. INTRODUCTION
The origin of dark matter is on the basis of the aim to elucidate the existence of new low-mass, weakly interacting particles from a theoretical, phenomenological and experimental point of view. The light particle that has gathered the most attention is the axion, predicted by Weinberg and Wilczek [1], and Wilczek [2] in the Peccei and Quinn mechanism [3] to explain the absence of parity and temporal invariance violations induced by the QCD vacuum. The axion is one of the more interesting candidates to make the dark matter of the universe, and the axion potential, that determines the dynamics of the axion field, plays a fundamental role in this context. The QCD axion model relates the topological susceptibility χ T with the axion mass m a and decay constant f a through the relation χ T = m 2 a f 2 a . The axion mass is, on the other hand, an essential ingredient in the calculation of the axion abundance in the Universe. Therefore, a precise computation of the topological properties of QCD and of their temperature dependence becomes of primordial interest in this context. Understanding the role of the θ parameter in QCD and its connection with the strong CP problem is one of the major challenges for high energy theorists [4].
The calculation of the topological susceptibility in QCD is already a challenge, but calculating the complete potential requires a strategy to deal with the so called sign problem, that is, the presence of a highly oscillating term in the path integral. In fact euclidean lattice gauge theory, our main non-perturbative tool for studying QCD from first principles, has not been able to help us much because of the imaginary contribution to the action coming from the θ term, that prevents the applicability of the importance sampling method [5]. This is the main reason why the only progress in the analysis of the finite * Corresponding author: eduroyo@unizar.es temperature θ dependence of the vacuum energy density in pure gauge QCD, outside of approximations, reduces to the computation of the first few coefficients in the expansion of the free energy density in powers of θ [6], and the situation in full QCD with dynamical fermions is, on the other hand, even worse [7][8][9][10][11][12].
Much experience has been developed in the last years by our group in this field, both in the elaboration of efficient algorithms to simulate systems with a theta-vacuum term overcoming the severe sign problem [13,14], as well as in the application of these approaches to the computation of the vacuum energy density and topological charge density for several interesting physical systems [15][16][17][18][19]. Our purpose is to take advantage of this experience to apply these approaches to the computation of the θ dependence of the QCD vacuum energy density.
As a first step in this ambitious program we present in this paper an analysis on the θ dependence of a toy model for QCD, the Schwinger model, on the lattice. Strictly speaking, the Schwinger model in the continuum is not asymptotically free, as QCD, since it is superrenormalizable and the Callan-Symanzik β-function vanishes. However, in the lattice version, since the continuum coupling is dimensionful, the continuum theory is reached at infinite inverse square gauge coupling β = 1/e 2 a 2 , much in the same way as four-dimensional asymptotically free gauge theories such as QCD. Furthermore the model is confining [20], exactly solvable at zero fermion mass, has non-trivial topology and shows explicitly the U A (1) axial anomaly [21] through a non-vanishing value of the chiral condensate in the chiral limit in the one-flavor case. These are basically the reasons why this model has been extensively used as a toy model for QCD.
For two dimensional systems such as the Schwinger model with a θ term, there exist numerical methods such as the Hamiltonian method [22,23] and the Grassmann tensor renormalization group method [24] that have been applied successfully. However, such methods are currently only applicable to two-dimensional systems, whereas our aim is to test a method that should, in principle, be applicable also to four-dimensional theories such as QCD.
The paper is organized as follows; in Section II we summarize some relevant features of the Schwinger model with topological term. Since our second proposal to analyze physical systems with a topological term in the action [14] has been found to be particularly well suited to bypass the sign problem in asymptotically free gauge theories, we decided to apply it, and Sec. III contains a brief review of the main steps of the method. In Sec. IV we give details on the lattice setup and the computer simulations. Sec. V shows our results for the topological charge density as a function of θ at several fermion masses and gauge couplings, and finally in Sec. VI we report our conclusions.

II. THE MASSIVE SCHWINGER MODEL WITH A θ TERM
The Schwinger model is Quantum Electrodynamics in 1+1-dimensions [25]. The euclidean continuum action reads where m is the fermion mass and e is the electric charge or gauge coupling, which has the same dimension as m.
After a simple rescaling of the fields the action can be written as where and γ µ are 2 × 2 matrices satisfying the algebra At the classical level this action is invariant in the chiral limit under the U A (1) global transformations leading to the conservation of the axial current However the axial symmetry is broken at the quantum level because of the axial anomaly. The divergence of the axial current is with ǫ µν the antisymmetric tensor, and therefore does not vanish. The axial anomaly induces a topological θ term in the action of the form where the topological charge Q = 1 4π d 2 xǫ µν F µν (x) is an integer.
The purpose of this paper is to analyze the θ dependence of the model described by the action (2)+(8) A simple analysis of this model on the lattice suggests that it should undergo a phase transition at some intermediate fermion mass m and θ = π, even at finite lattice spacing. Indeed the lattice model is analytically solvable in the infinite fermion mass limit (pure gauge twodimensional electrodynamics with topological term) [26], and it is well known that the density of topological charge approaches a non-vanishing vacuum expectation value at θ = π for any value of the inverse square gauge coupling β, exhibiting spontaneous symmetry breaking. On the other hand by expanding the vacuum energy density in powers of m, treating the fermion mass as a perturbation [27], one gets for the vacuum expectation value of the density of topological charge the following θ dependence: with Σ the vacuum expectation value of the chiral condensate in the chiral limit and at θ = 0 (Σ = e γe e/2π 3/2 in the continuum limit), and χ P and χ S the pseudoscalar and scalar susceptibilities respectively. Equation (10) shows how the Z 2 symmetry at θ = π is realized order by order in the perturbative expansion of the topological charge in powers of the fermion mass m, and therefore a critical point separating the large and small fermion mass phases is expected. Indeed the model was analyzed in the continuum by Coleman in [28], where he conjectured the existence of a phase transition at θ = π, and some intermediate fermion mass m separating a "weak coupling" phase ( e m << 1), where the Z 2 symmetry of the model at θ = π is spontaneously broken, from a "strong coupling" phase ( e m >> 1) where the Z 2 symmetry is realized in the vacuum. This conjecture was corroborated in [22,23] using the lattice Hamiltonian approach with staggered fermions, and more recently in [24] using the Grassmann tensor renormalization group and Wilson fermions.

III. COMPUTING THE ORDER PARAMETER AS A FUNCTION OF θ
To compute the θ dependence of the density of topological charge we use the approach proposed in reference [14]. The only assumption in this approach is the absence of phase transitions at real values of θ except at most at θ = π. The method is based in extrapolating a suitably defined function to the origin. This function turns out to be very smooth in all the cases considered up to now [15][16][17][18], and this makes us confident on the whole procedure. Here we summarize the main steps.
From numerical simulations of our physical system at imaginary values of θ = −ih (real values of h), which are free from the severe sign problem, we compute the density of topological charge q(−ih) as a function of h, and introduce the following functions: The procedure to find out the density of topological charge at real values of θ relies on scaling transformations [14]. We define the function y λ (z) as For negative values of λ, the function y λ (z) allows us to calculate the order parameter tanh h 2 y (z) below the threshold z = 1. If y (z) is non-vanishing for any positive z, 1 then we can plot y λ /y against y. Furthermore, in the case that y λ /y is a smooth function of y close to the origin, then we can rely on a simple extrapolation to y = 0. Of course, a smooth behavior of y λ /y cannot be taken for granted; however no violations of this rule have been found in the exactly solvable models.
The behavior of the model at θ = π can be ascertained from this extrapolation. At θ = π the model has the same Z 2 symmetry than at θ = 0. We can define an effective exponent γ λ by As z → 0, the order parameter tan θ 2 y cos θ 2 behaves as (π − θ) γ λ −1 . Therefore, a value of γ λ = 1 implies spontaneous symmetry breaking at θ = π. A value between 1 < γ λ < 2 signals a second order phase transition, and the corresponding susceptibility diverges. Finally, if γ λ = 2, the symmetry is realized (at least for the selected order parameter), there is no phase transition and the free energy is analytic at θ = π. 2 We can take the information contained in the quotient y λ y (y), and calculate the order parameter for any value of θ through an iterative procedure [14]. The outline of the procedure is the following: i. Beginning from a point y (z i ) = y i , we find the value y i+1 such that y λ = y i . By definition, y i+1 = y e −λ 2 z i .
ii. Replace y i by y i+1 , and start again.
The procedure is repeated until enough values of y are know for z < 1 (see Fig. 1). This method can be used for any model, as long as our assumptions of smoothness and absence of singular behavior are verified during the numerical computations. The reliability of our approach in practical applications is better when the following conditions are met: 2 Other possibilities are allowed, for instance, any γ λ > 1, γ λ ∈ N leads to symmetry realization for the order parameter at θ = π and to an analytic free energy. If γ λ lies between two natural numbers, p < γ λ < q, p, q ∈ N, then a transition of order q takes place.
a. y(z) takes small values for values of z of order 1.
b. The dependence on y of the functions y λ /y and γ λ is soft enough to allow a reliable extrapolation.
In the one-dimensional Ising model within an imaginary magnetic field these two properties are realized in the low temperature regime [14], and the two and threedimensional models also show a very good behavior in this regime [17]. Indeed the relevant feature, at least in what concerns point a, is that, at low temperatures, the magnetic susceptibility at small values of the real external magnetic field takes small values. In the more interesting case of asymptotically free models, the analogue of the magnetic susceptibility is the topological susceptibility, and it is well known that topological structures are strongly suppressed near the continuum limit. Therefore, and on qualitative grounds, we expect a good implementation of our method in the Schwinger model or in QCD, at least close enough to the continuum limit. In fact, concerning asymptotically free models, the method was successfully applied to the analysis of the continuum θ−dependence of CP 9 [15], showing a very good realization of conditions a and b. In the more general cases we should find out whether the model agrees with these two conditions or not, and pleasant surprises are not excluded [16], [18].

IV. DETAILS OF THE SIMULATION
We use the lattice version of the continuum action (9) with staggered fermions and standard Wilson form for the pure gauge part. It reads as follows, where the notation is standard. The compact gauge variable U µ (n) is related to the non-compact gauge field A µ (n) in the usual way with a the lattice spacing, and the local topological charge q(n) is given by An important point is that this charge is quantized, and therefore the partition function of the model has exact 2π periodicity in θ, as is the case in the continuum theory. 3 We will analyze in what follows the model given by action (15), taking the square root of the fermion determinant in order to describe only one flavor. There is ample evidence that this procedure leads to the correct physics in the continuum limit, including the effects of the anomaly. For example, the Microcanonical Fermion Average (MFA) approach [30] was applied years ago to simulate the one-flavor Schwinger model on the lattice at θ = 0, using the standard Wilson action for the gauge field and staggered fermions. The results [31] reproduce the exact value of the chiral condensate in the chiral continuum limit up to 3 decimal places. 4 As explained in Sec. III, in order to find the dependence on θ of the density of topological charge q, we need to compute the expected value of q for imaginary values of θ. In this case standard Monte Carlo algorithms work well, and we can sample the distribution e −S generated by (15) with any of these methods. We have used a standard Metropolis approach, trying to update each link sequentially in every sweep.
We want to use an exact Monte Carlo method, and the fermionic part of the action forces us to recompute the whole fermion determinant at each attempt to update one link, that is, N times each sweep. Indeed, this is the most expensive part of the algorithm. All our lattices are of size 16 × 16, and we computed the eigenvalues of the fermion matrix with the GNU Scientific Library, taking advantage of the standard even-odd decomposition of the staggered fermions Dirac operator. The simulations have been run at the U-LITE computer facility at the INFN National Laboratories of Gran Sasso.
We present in the following section results for several masses m, gauge couplings β, and fields h (imaginary θ). At each point of the parameter space, we run the algorithm and take up to 100k measurements, each one made every ten sweeps. Between 5-10% of the initial configurations are discarded for thermalization. The errors of the expected values are estimated by a standard jackknife binning.

V. RESULTS
Our results for the exponent γ are summarized in Figures 2 and 3. As is apparent in Fig. 2, the behavior at fixed β is very different as we vary the fermion mass. The data corresponding to m = 0.5 lie essentially on top of the analytic m = ∞ curve (that is, pure gauge theory [26]), and therefore this mass corresponds to the phase with broken symmetry at θ = π. On the other hand, the data for both the m = 0.05 and the m = 0 case extrapolate to a value clearly above 1, indicating symmetry restoration at θ = π, although our data are not precise enough to make a definite statement on the value of γ. In Fig. 3 we present the results at fixed m = 0 for the various coupling constants we have studied. We see as before a clear extrapolation to a value of γ above 1, 5 in stark contrast with the pure gauge theory case. We have also been able to extract the full dependence of the order parameter as a function of the angle, q(θ). In Fig. 4 we show details of the fit y versus y λ for a particular value of the parameters, in order to give an idea of the precision of our data. This will be followed by the iterative procedure depicted in Sec. III in order to produce the curve q(θ). Regarding the error estimation, it is impossible to follow the error propagation from the values of q at imaginary θ until the final result of q(θ). To overcome this impasse we proceed in the following way: first we generate 20 sets of fake data for q(h) having the same mean value and distribution of the actual Monte Carlo data; then we compute, following the same scheme (fit of y vs. y λ and iterative procedure), 20 realization of q(θ); from the distribution of these values around the curve computed from the real data, we can infer the error associated to each value of q(θ), which will be shown in the following figures as a shaded band. In Fig. 5 we present q(θ) at β = 3 for two masses in the symmetry restored phase, as well as at β = 2 and In Fig. 6 we show the results for m = 0 and the three different values of the coupling constant we have simulated. We can clearly see the restoration of the symmetry as we approach θ = π. In Fig. 7 we show, for β = 3.0 and m = 0, the order parameter q(θ) in the vicinity of θ = π. Fitting q(θ) near θ = π in the symmetry restored phase allows us to extract the exponent (π − θ) ǫ , which is related to γ by ǫ = γ − 1. 6 We present in Table I our results for ǫ.
To finish this Sec. we want to discuss a little bit more on the results for the massless Schwinger model reported in Fig. 6. It is well known that the continuum formulation of the massless Schwinger model shows no θ dependence, because the θ term in the action can be canceled by an anomalous chiral transformation which does not change the fermion-gauge action if the fermion mass vanishes. Hence the non-trivial θ dependence of the density of topological charge shown in Fig. 6 may seem surprising. However, the massless staggered Dirac operator does not have exact zero-modes, and therefore, for a given gauge configuration, a nonzero value of the quantized topological charge Q does not imply the existence of a corresponding number of zero-modes in the staggered Dirac operator, as would be the case, for example, with the overlap Dirac operator. What we should expect instead is that, as we approach the continuum limit, the topological charge density vanishes. This is indeed what seems to happen, as is suggested by the results of Fig. 8.

VI. CONCLUSIONS AND OUTLOOK
All our results are compatible with the standard lore on this model, and in particular with Coleman's conjecture on the existence of two distinct phases at θ = π, a symmetry breaking phase at large mass, and a symmetry restored phase at small mass.
Our simulations are a proof of concept, and are not extensive enough to determine precisely the position of the critical mass at θ = π or its properties in detail. But the important point is that we have succeeded in calculating the full dependence of the order parameter in θ in a gauge theory with fermions and a quantized topological charge, using a method that should, in principle, work also in higher dimensional theories.
The aim of this calculation was to test the method developed in [14] in a fermionic gauge theory, as a first step towards applying it to QCD with a θ term. In light of the excellent results obtained, we expect the method to be applicable also in this case.