Benchmark results in the 2D lattice Thirring model with a chemical potential

We study the two dimensional lattice Thirring model in the presence of a fermion chemical po- tential. Our model is asymptotically free, contains massive fermions that mimic a baryon and light bosons that mimic pions. Hence it is a useful toy model for QCD, especially since it too suffers from a sign problem in the auxiliary field formulation in the presence of a fermion chemical potential. In this work we formulate the model in both the world line and fermion-bag representations and show that the sign problem can be completely eliminated with open boundary conditions when the fermions are massless. Hence we are able accurately compute a variety of interesting quantities in the model, and these results could provide benchmarks for other methods that are being developed to solve the sign problem in QCD.


I. INTRODUCTION
Traditional lattice calculations of quantum field theories often encounter sign problems in the presence of a chemical potential. An excellent example is Quantum Chromodynamics (QCD), where it is impossible to accurately compute quantities at a non-zero baryon density especially at low temperatures [1]. Over the past decade, ideas like the complex Langevin approach [2] and the Leftchetz thimble approach [3] have been proposed as potential solutions to sign problems including QCD. When these methods are tested on simple models where exact results are available [4][5][6], we not only find potential pitfalls of the methods but also learn new directions to avoid them [7][8][9]. While these ideas have also been able to capture some of the qualitative features of more complex field theories [10,11], in these cases the numerical results are not always compared with benchmark calculations obtained with other methods where the errors can be controlled. An exception to this has been studies of bosonic field theories at finite densities where a controlled Monte Carlo algorithm in the world line representation free of sign problems is available [12,13]. Producing such benchmark calculations that truly test the method, especially in fermionic quantum field theories with a sign problem and similar to QCD in other aspects, would be helpful and is the main motivation behind our work.
Recently the Lefshetz thimble program got a boost when it was shown that it may be possible to use holomorphic flow in complex field space to sample multiple thimbles rather than perform calculations on a single thimble as was done in the past [14]. The focus has also turned to lattice Thirring models as a prototype example of the physics of QCD [14,15]. This model has also been studied earlier in higher dimensions using stochastic quantization [16]. Also, the recent work has computed the average fermion number N on small but fixed spatial size L X as a function of the chemical potential, which is much more sensitive to the important physical scales in the problem, as compared to local densities on large space-time lattices. As shown schematically in Fig. 1, at low temperatures (or large L t ) the plot of N as a function of the chemical potential µ is expected to show a series of jumps at critical values of the chemical potential, say µ 1 , µ 2 , ... where the average particle number jumps to N 1 , N 2 , .... The values of µ i 's and N i 's are related to the physical scales of the problem like binding energies and scattering lengths and should become harder to compute due to sign problems, especially when µ i and N i become large. Encouraged by the fact that some of these quantitative features may be within reach, recently the efforts have turned towards speeding up the calculations on larger lattices using machine learning algorithms [17]. It would indeed be exciting if this program is successful.
The motivation for our work is to help this program by accurately computing the µ i 's and N i 's for a specific two dimensional lattice Thirring model constructed with staggered fermions. Our model is asymptotically free and a continuum limit can be defined at zero cou- FIG. 1: A schematic plot of the particle number as a function of chemical potential for a fixed spatial size. We propose that the values of µi and the corresponding Ni's are easily calculable and can be used as benchmark quantities to validate a method that claims to solve a sign problem.  fermion in the theory is massive and mimics a baryon,  while bosonic excitations made with fermion-antifermion  pairs are massless and mimic pions. Thus, the similarities of our model with QCD is striking. Of course the  ground state does not break any symmetries and the pions are not really Goldstone bosons as was explained by Witten long ago [18], but the fermion mass generation is dynamical like in QCD and from the point of view of sign problems the bosons being lighter than the fermion is also similar to QCD. Interestingly, we can solve the model both in the fermion world line method and the fermion bag approach. In the world line approach we argue in this work that the sign problem is absent with open boundary conditions and zero fermion mass. Thus, in this limit we are able to study large lattices and can accurately compute the the critical µ i 's and N i 's. These could provide helpful benchmark to test new ideas that claim to solve sign problems in problems similar to QCD. Our paper is organized as follows. In section II we discuss the model we study and the various types of representations that can be used to solve it. In particular we show why the model in the massless limit with open boundary conditions has no sign problem in the worldline formulation. In section III we discuss our Monte Carlo methods, especially the worm algorithm to update the worldline representation and the fermion bag algorithm. In section IV we discuss the results we have obtained. In particular we define the observables we measure and discuss our results in a variety of parameter ranges. We present our conclusions in section V.

II. THE MODEL
The lattice action of the model we study is given by (2.1) where the the matrix M is the massless staggered fermion matrix defined as where µ is the chemical potential, m is the fermion mass and η x,µ are the usual staggered phase factors (η x,0 = 1 and η x,1 = (−1) x1 ). The four-fermion coupling U can be interpreted as a current-current interaction on neighboring sites and hence the name "lattice Thirring model". When m = 0 the model contains the well known U (1) chiral symmetry of staggered fermions. In the discussion below L X denotes the number of spatial sites and L T denotes the number of temporal sites in our two dimensional square lattice. Further, we always use anti-periodic boundary conditions in time, but study the effects of periodic, anti-periodic and open boundary conditions in space.
This model has a long history and has been studied extensively in three space-time dimensions in the auxiliary field formulation [19,20] and the fermion bag approach [21,22]. In three dimensions the model with m = 0 has two phases: a weak coupling phase with massless fermions and a strong coupling phase with spontaneously broken U (1) chiral symmetry, massive fermions and light pions. These phases are separated by a second order critical point, whose properties were studied in the earlier work. In two dimensions this critical point moves to the origin and the massless weak coupling phase disappears. Further, since a continuous chiral symmetry cannot break in two dimensions, the massive fermion phase becomes critical. Thus, the two dimensional model contains massive fermions and critical bosons, where the mass of the fermion can be used to set the lattice spacing. The continuum limit is taken by tuning U towards the origin. As far as we know these features of the two dimensional model with m = 0 were never studied using the Monte Carlo method even at µ = 0 where there is no sign problem. The similarity of the model with QCD makes it an interesting toy model for studies at non-zero chemical potential. At a large value of m, this was done recently in two space-time dimensions [23].

A. The Auxiliary Field Representation
The traditional approach to solve these models is by rewriting the partition function using an auxiliary field formulation so that it can be tackled by the Hybrid Monte Carlo algorithm. More explicitly where in the last step we have introduced a compact auxiliary field 0 ≤ A x,ν < 2π associated with the bonds of the lattice and the auxiliary field action 4) is now a Gaussian in the Grassmann fields. The Dirac matrixM x,y is defined as ν η x,ν 2 e iAx,ν +µδν,0 δ x+ν,y − e −iAy,ν −µδν,0 δ x,y+ν (2.5) and the parameters U and m are related to g and m through the relations (2.6) Here I 0 and I 1 are the Bessel function and the first modified Bessel function. The sign problem in the auxiliary field representation can be traced to the fact that Det(M + m ) does not have any symmetries and can be complex when µ = 0, like in QCD.

B. The Fermion Bag Representation
Can ideas of fermion bags help solve the sign problem present in the auxiliary field approach? In this approach we do not introduce the usual auxiliary fields, but try to regroup fermion worldlines differently. Unfortunately, this regrouping is not unique and needs some thought. One possible regrouping introduced earlier for µ = 0 case is based on introducing a new set of variables, the dimers d x,ν for nearest neighbor interactions and monomers n x for the mass terms [21]. This naturally emerges when we expand the Grassmann exponential of the mass and interaction terms: We then interpret the expression As an illustration we show a possible configuration of dimers and monomers on a 4 × 4 block of lattice sites in Fig. 2  When µ = 0, since the matrices W ([f ]) are always anti-symmetric, Det(W ([f ]) ≥ 0 and the sign problem is solved. However, in the case of µ = 0 this property no longer holds and the determinants can be negative. This may seem surprising since in two space-time dimensions the fermion permutation sign is absent due to the fact that fermions cannot cross each other. In our model fermions have a flavor and they can change flavors while hopping. This is encoded in the staggered phase factors and this leads to a sign problem. Empirically we discovered that this remaining sign problem depends on the boundary conditions. While the sign problem is present with both periodic and anti-periodic boundary conditions, it is absent with open boundary conditions. This also means that on large space-time lattices with L X = L T the sign problem essentially disappears, but for asymmetric lattices it can reemerge. In the most interesting case for our studies, where we fix the spatial lattice size L X , and study very large values of L T the sign problem can become severe with periodic and antiperiodic boundary conditions.

C. World Line Representation
In order to get a better understanding of the origin of the sign problem in our model, we look at the representation of the fermion determinant Det(W [f ]) inside free fermion bags as a sum over their world lines. This representation can be found by expanding the determinant back into the Grassmann integral form, This product can be represented in terms of directed fermion link variables l x,±ν = 0, ±1, where +1 represents the term χ x χ x±ν and −1 represents the term χ x+ν χ x . The determinant is replaced with a sum over all configurations of directed links.
Configurations of links only have a nonzero weight when oneχ and one χ are chosen at each site. Thus each site must have one directed link pointing into it and one pointing out of it. The links will therefore form closed loops.. In Fig. 3 we show two valid configurations with the directed links represented as arrows pointing from χ to χ. The weight of a configuration of fermion world lines is given by the product of the weights in Eq. 2.10 and a factor −1 for every closed loop arising from a reordering of χ x andχ x to match the ordering of the measure.
where N loops is the number of closed loops formed by the directed links. It is easy to verify that there are valid configurations with a negative weight. For example, the configuration on the left in Fig. 3 has a positive weight, but the configuration on the right has a negative weight.
Let us now prove that the sign problem disappears with open boundary conditions in the massless limit because configurations with a negative sign are absent at the worldline level. The weight of a configuration can be written as the product of the weights of the closed loops of fermion links It is therefore sufficient to show that all loops that can exist in a configuration have positive weight. Consider first a loop that does not wrap around the volume. Note that by starting from the trivial loop that visits two neighboring sites, we can construct any nonwrapping loop using the two deformations depicted in  The second deformation inverts a corner and introduces a single site inside the loop. This does change the sign of the loop. Thus, any non-wrapping loop can be negative only if it encloses an odd number of sites. But in the massless limit a configuration with such a loop will have zero weight, since monomers are not allowed, dimers always take away two sites and all free fermion loops will touch an even number of sites. Given that the trivial loop has a positive sign, any allowed non-wrapping loop will have a positive sign. Similarly, any loop enclosing an odd number of sites has a negative sign.
Thus, any negative sign fermion loops in the massless theory must arise through loops that wrap around the temporal direction. Note that with open boundary conditions spatial winding is also forbidden. We can again construct any temporal wrapping loop by starting from a loop that goes straight in time without hops and deforming it using the two deformations discussed above. This time a negative sign in the loop can be introduced if an odd number of sites cross the fermion line during this deformation. An example of such a negative signed loop is shown in Fig. 5  the boundary. However, with open boundary conditions such temporal loops will create regions on the left and right with an odd number of sites. This is forbidden in the massless limit for the same reasons as outlined above for non-wrapping loops.
With periodic and anti-periodic boundary conditions we can have other more complicated loops as shown in Fig. 5 on the right. Thus, the sign problem can be completely eliminated by open boundary conditions in the spatial direction. This feature of the world line formulation is well known and specific to two dimensional models [24][25][26]. In higher dimensions the argument for the positivity of all fermion worldline configurations fails and significant cancellations between world line configurations will be necessary for alleviating the sign problem. The fermion bag approach can be helpful in this regard [23].

III. MONTE CARLO UPDATES
Monte Carlo methods for updating both the worldline representation and the fermion bag representations are by now well developed [27][28][29][30][31][32]. We use a worm algorithm to update the fermion lines and dimers, using updates like the one illustrated in Fig. 6. To begin an update, we suggest randomly changing some fermion link l x,ν . The update is then accepted with the absolute value of the weight given in Eq. (2.11). If the link is changed two defects are generated in the lattice configuration which are allowed. The defects are the head and tail of the worm. The head of the worm then propagates by updating the neighboring links. When the head returns to its tail, the worm closes, the defects disappear and the update is complete. The various steps of how the defect propagates are shown in shown Fig. 6. The configuration of dimers may also be updated during the worm update. When this is done we have to use the weights of including or removing a dimer. Fig. 7 shows the steps for an update that changes the dimer number.
In contrast to the worm algorithm, we sample fermion bag configurations using a local Monte Carlo update that involves adding or removing dimers or pairs of monomers. Each proposal is accepted with the probability where the new configuration is denoted with primed variables. The fact that one has to use ratios of fermion determinants that are non-local helps in reducing autocorrelation times. We also can update large regions of space-time by using a background field method used recently in [32]. The sampling is made more efficient with a move that switches the places of a monomer and a dimer if the two are on neighboring sites. Since the weights of the two configurations are the same this update is very quick.

IV. NUMERICAL RESULTS
In this work we compute three observables in order to understand the physics of our model. The first is the chiral condensate susceptibility χ, defined by the relation We can use it to understand the physics of bosonic excitations in our model. We also compute the chiral charge  winding number susceptibility, defined by the relation where S and S are surfaces orthogonal to the direction α.
In the thermodynamic limit, the winding number susceptibility helps us understand the status of chiral symmetry as we explain below. In the world line representation the chiral charge can be defined by the relation which means the susceptibility is simply since the chiral charge is conserved on each configuration. Finally we measure the average fermion number using the relation N f = x∈S J 0,x where the fermion number current is given by and S is a surface perpendicular tot. In the worldline representation again the fermion number is straight forward to calculate and is given by In our definition the fermion number is normalized to count both the Dirac and flavor degrees of freedom from a continuum limit perspective. Using these observables we first focus on the physics of our model at µ = 0 in order to bring out the similarities to QCD. As we mentioned earlier, unlike in QCD the U (1) chiral symmetry of the model cannot break in two dimensions. However, the lightest boson in the model is critical (i.e., it is massless but is not a Goldstone boson). Hence when L X = L T = L we expect the chiral condensate susceptibility to scale as χ ∼ L 2−η for large values of L. The exponent η depends on U like in the usual critical phase of the two dimensional XY model. At infinite U since the Thirring model becomes a closed packed dimer model and we expect η = 0.5 [33]. When U = 0 the susceptibility diverges logarithmically with L and hence η = 2. Our results reproduce this and show how the exponent changes continuously between these two limits. In table I  expected to vanish because the chiral charge cannot wind across the spatial boundaries. However, when the phase is critical like in our model it is expected to be go to a constant in the thermodynamic limit. Our results are consistent with this expectation. The values we measured for Q 2 χ at L = 256 are given in table I. These values are found using open boundary conditions Further we find that Q 2 χ = 0.25 at U = 0, and grows monotonically to the value of roughly 1.2 at U = ∞. All this is consistent with the fact that the bosonic sector of our theory is critical.
In contrast to the bosons, fermions are massive for all values of U > 0. We compute the fermion mass m f as a function of U using large square lattices (L X = L T = L) as follows. In the thermodynamic limit we expect the average fermion density n = N f /L X to be zero when µ ≤ m f and rise linearly according to the relation for µ ≥ m f . This behavior should also be an excellent approximation for sufficiently large lattices. To demonstrate this we show our results for n at U = 0.3 with open boundary conditions in Fig. 8. Selected values of this data are also tabulated in table II as a benchmark for future calculations. As we can see for L = 10 the curve does not show the expected non-analyticity, but for L = 40 and L = 64 the curves show it clearly. We can fit our data to the linear form which is shown as the dashed line in the fit. In table I we report the value of m f found using this method for several values of U . We used lattices of size of L = 64, except for U = 0.1, where the lattice size used was L = 128. The dynamical generation of fermion mass is an interesting feature of our model. While similar to the phenomenon of chiral symmetry breaking in QCD, the actual dynamical breaking of continuous symmetries is forbidden in two dimensional models. Nevertheless a fermion mass can be generated and a massless boson with critical correlations can arise [18]. Finally we note that fourfermion couplings are expected to be marginal in two dimensions and in our case it also happens to be marginally relevant (i.e., asymptotically free). Thus, at small U the fermion mass m f is expected to vanish according to the relation where b 0 = 16 is the one-loop coefficient of the β function. Figure 9 shows the fermion mass values and compares it against the expected behavior. For purposes of illustration we use C = 0.49. With these small masses lattice volumes up to V = 1024 × 1024 were necessary. It is well known that such asymptotic scaling fits don't work very well unless very large lattices are used [34]. Here we just use it to illustrate that the the fermion mass qualitatively does become exponentially small as U becomes small. Next we turn to the physics of finite chemical potential. We first consider a small spatial lattice of L X = 6 and L T = 48 and study the sign problem in the traditional auxiliary field approach with periodic boundary conditions and compare it with the sign problem in the fermion bag approach with both periodic and antiperiodic boundary conditions. In Fig. 10 we plot the average sign as a function of the chemical potential in the auxiliary field approach at U = 0.3 (left) and compare it with that of the fermion bag approach (right). We first wish to learn where the sign problem becomes severe. In the auxiliary field approach the sign becomes severe around µ ≈ 0.4, while in the fermion bag approach with anti-periodic boundary conditions it becomes severe around µ ≈ 0.55. In the fermion bag approach with periodic boundary conditions the sign problem is never severe, although it is enhanced both at µ ≈ 0.4 and then again at µ ≈ 0.9. Can we correlate this behavior with some physics?
Let us now explore how the fermion chemical potential "dopes" the system with fermions. We again focus first on a small lattice, L X = 6, L T = 48 at U = 0.3. In table III we present all of our results for the total fermion number as a function of the chemical potential for open(left), anti-periodic(center) and periodic boundary conditions(right). In Fig. 11 we plot these results  along with the results for free fermions as solid lines. Due to the flavor degeneracy of staggered fermions we expect all states to be at least doubly degenerate. With open periodic boundary conditions this means all jumps must be in steps of two. This is what is observed. With periodic and anti-periodic boundary conditions there is a symmetry between left and right moving particles. With periodic boundary conditions a zero momentum state is allowed which is non-degenerate, hence the first jump in N near µ ≈ 0.4 is only by two. However, the second jump near µ ≈ 0.9 is by a factor of four since now nonzero momentum states are excited and each state is doubly degenerate due to the two fermion flavors. With antiperiod boundary conditions the lowest energy state already has momentum and hence again should have fourfold degeneracy. This is clearly seen as a jump of four in the free theory around µ ≈ 0.5. Surprisingly, in the interacting theory this degeneracy of the lowest energy state seems to be broken. We attribute this to the fact that bound state bosons with zero momentum can emerge. The next momentum state is non-degenerate for L X = 6 since effectively the lattice size is halved for staggered fermions. This remains unchanged for the interacting theory as well and two additional states are added when µ > 1.
Note that the first step to N f = 2 for both open and periodic boundary conditions occurs around µ ≈ 0.4. This coincides with the point where the sign problem becomes severe in the auxiliary field approach, and is somewhat enhanced in the fermion bag approach. The sign problem in the fermion bag approach disappears for large values of µ until around µ ≈ 0.9 where there is the second jump of four in the periodic case. The sign problem in the auxiliary field approach on the other hand never recov-   Fig. 13. ers. In the case of anti-periodic boundary conditions the severity of the sign problem coincides with the additional plateau at N f = 2 which is absent in the free theory as discussed above. While these correlations between sign problems and the underlying physics are not surprising, the fact that energies and degeneracies of the lowest lying states can be influenced by boundary conditions and interactions on small lattices offers an excellent opportunity for methods that claim to solve the sign problems to reproduce them.
Since the sign problem is absent with open boundary conditions we can use it to study the behavior of N f   Fig. 13. on large asymmetric lattices (L X = L T ) so as to understand the physics of fermion doping at a fixed L X . One of the main results that our model shares with QCD is that fermions become massive entirely due to interaction effects and the value of the chemical potential where the first jump in N f occurs will be this finite size fermion mass m L X f . In order to see the effects of interactions we plot N f as a function of µ in the free theory ( Fig. 12) and in the interacting theory with U = 0.3 ( Fig. 13) both with open boundary conditions. Selected data points have also been tabulated in Tables. IV,V and VI for benchmark purposes.
We study three different lattice sizes L X = 12 (left) L X = 16 (center) and L X = 32 (right). For each of these lattices we study the effects of increasing L T . Note that the critical value of µ where the first jump to N f = 2 occurs, shifts to lower values as L X increases in the free theory. We expect this value to vanish in the large L X limit since fermions are massless. However, in the interacting theory we note the jump change in the critical value is smaller and should approach 0.183(1) (see table I) as L x becomes large. Also the jump becomes sharper as the anisotropy (value of L T ) is increased and approaches a step function as expected. To quantify the value of m L X we can also extract the finite size boson mass m L X b . These values are also given in table VII for L X = 12 and 32. We find that while m L X f increases sharply with U , m L X b decreases mildly.

V. CONCLUSIONS
In this work we have studied the 1 + 1 dimensional lattice Thirring model with staggered fermions at both zero and finite densities. We showed that the model is free of sign problems in the massless limit when open boundary conditions are used. In this case we used the worldline formulation to study the model. In the case of periodic and anti-periodic spatial boundary conditions the sign problem is mild on square lattices but becomes severe when on asymmetric lattices. However, the fermion bag formulation seems to alleviate the problem except at critical values of the chemical potential where fermion number jumps. We provide accurate estimates for the total particle number as a function of the chemical potential for a few lattice sizes. Our results could be used as a benchmark for future studies by other methods that attempt to solve the sign problem.