Absence of chiral symmetry breaking in Thirring models in 1+2 dimensions

The Thirring model is an interacting fermion theory with current-current interaction. The model in $1+2$ dimensions has applications in condensed-matter physics to describe the electronic excitations of Dirac materials. Earlier investigations with Schwinger-Dyson equations, the functional renormalization group and lattice simulations with staggered fermions suggest that a critical number of (reducible) flavors $N^{\mathrm{c}}$ exists, below which chiral symmetry can be broken spontaneously. Values for $N^{\mathrm{c}}$ found in the literature vary between $2$ and $7$. Recent lattice studies with chirally invariant SLAC fermions have indicated that chiral symmetry is unbroken for all integer flavor numbers [Wellegehausen et al., 2017]. An independent simulation based on domain wall fermions seems to favor a critical flavor-number that satisfies $1<N^{\mathrm{c}}<2$ [Hands, 2018]. However, in the latter simulations difficulties in reaching the massless limit in the broken phase (at strong coupling and after the $L_s\to\infty$ limit has been taken) are encountered. To find an accurate value $N^{\mathrm{c}}$ we study the Thirring model (by using an analytic continuation of the parity even theory to arbitrary real $N$) for $N$ between $0.5$ and $1.1$. We investigate the chiral condensate, the spectral density of the Dirac operator, the spectrum of (would-be) Goldstone bosons and the variation of the filling-factor and conclude that the critical flavor number is $N^{\mathrm{c}}=0.80(4)$. Thus we see no chiral symmetry breaking in all Thirring models with $1$ or more flavors of ($4$-component) fermions. Besides the artifact transition to the unphysical lattice artifact phase we find strong evidence for a hitherto unknown phase transition that exists for $N>N^{\mathrm{c}}$ and should answer the question of where to construct a continuum limit.


I. INTRODUCTION
The Thirring model [4] in 2 space-time dimensions is integrable and in the massless limit even soluble [5,6]. The model in 3 space-time dimensions is of interest for various reasons, e.g. its close relationship to QED 3 [7][8][9][10] or its relevance in solid state physics, where it describes low-energy electronic properties of materials like graphene [11,12] or high-temperature superconductors [13,14]. In 3 dimensions the model is perturbatively non-renormalizable but can be renormalized in the limit of large flavor numbers N [7,[15][16][17]. Thus it provides a simple realization [18] of the concept of asymptotic safety [19]. In the large-N limit one finds an unbroken U(2N ) symmetry for every coupling strength. On the other hand, in the limit N = 1 /2 the Thirring model is equivalent to the Gross-Neveu model. The latter exhibits (for all N ) a second order phase transition from a symmetric gapless (massless) phase at weak coupling to a spontaneously broken gapped (massive) phase at sufficiently strong couplings 1 . We conclude that the Thirring model exhibits no chiral phase transition for large N but shows a second order phase transition at N = 1 /2. The question about the critical flavor number N c below which the Thirring model shows a chiral phase transition has been intensively discussed in the past. While early results obtained with functional methods or staggered lattice fermions range from N c = 2 to N c = ∞ * julian.johannes.lenz@uni-jena.de † bjoern.wellegehausen@uni-jena.de ‡ wipf@tpi.uni-jena.de 1 More precisely, the Thirring model with 1 irreducible 2-component Fermion flavor is the same as the Gross-Neveu model with 1 irreducible Fermion flavor. [8,[20][21][22][23][24][25][26][27][28][29][30], more recent lattice studies with chiral fermions favor smaller values of N c . In particular, based on simulations with massless (chiral) fermions we argued that the U(2N )-symmetry is unbroken for all integer flavor numbers N 1 [1]. For N = 1 the effective potential for the chiral condensate is almost flat at the origin such that we could not completely rule out the possibility that there is SSB for N = 1. In an interesting recent work Simon Hands applied domain-wall fermions to study the chiral condensate and masses of the (would-be) Goldstone bosons [3]. The results support our finding that N c is smaller than hitherto believed with the notable difference that he interprets his data as an evidence for 1 < N c < 2. In a recent explorative functional renormalization group (FRG) study with momentum-dependent couplings and the pseudo-spectral method the critical behavior of fourfermion theories [31] has been reconsidered. While a precise estimate for N c remains difficult in these elaborate FRG-studies, the new results are compatible with the lattice studies based on chiral fermions.
This work aims to solve the discrepancy between the results obtained with domain-wall and SLAC fermions. For that purpose we first performed simulations for 38 different non-integer values of N between 0.5 and 1.1 and calculated the corresponding chiral condensates. This way we already find strong evidence for a critical flavor number significantly lower than 1.0. However, due to the computational cost of the algorithm a reliable extrapolation to infinite volume is difficult. But with the help of a careful study of the (would-be) Goldstone spectrum and the spectral density of the Dirac operator we could not only assure unbroken symmetry at N ≥ 1.0 but also verify the proposed SSB U(2) → U(1) ⊗ U(1). We conclude that indeed there is a critical flavor number N c ≈ 0.80 which is considerably smaller than 1.0. However, the similarly accurate value for N c comes from studying the susceptibility of the four Fermi term in the Lagrangian which signals -besides the well-known transition to the artificial lattice phase -a new physical phase transition for all models with flavor numbers N ≥ N t = 0.78 (4). We argue that N t should be identified with N c . There is evidence that the new transition is of second order and can be used to construct a continuum limit of the lattice Thirring models. Interestingly, this new transition seems unrelated to any change of symmetry.
To summarize: All results of our simulations with SLAC fermions consistently show that chiral symmetry is not broken in all massless Thirring models with N = 1, 2, 3, . . . four-component fermions.
The paper is organized as follows: In the first section we recall relevant properties of the (reducible) Thirring model. For more details we refer to our earlier and much more detailed work [1], in which we investigated Thirring models with irreducible 2-component spinor-fields and with reducible 4 -component fields. In the present work we focus on the reducible and parity-even case considered in other works on the Thirring model in 3 dimensions. In the next two sections we present our lattice results for the chiral condensate and the spectral density -from which we extract a first estimate of the critical flavor number. Then we discuss the correlation matrix for interpolating operators for the scalar and pseudoscalar mesons. The simulation results for the meson spectra support the proposed symmetry breaking pattern of chiral symmetry. In the following section we present our simulation results for the expectation value of the interaction term ∝ (ψΓ µ ψ) 2 and the corresponding susceptibility. The expectation value is related to the mean filling factor of the fermions.
In appendix A we prove some useful properties of the spectral density and fermion Green function which follow from parity invariance of the reducible theory. Appendices B and C contain some technical details concerning numerical differentiation and our simulations.

II. THE THIRRING MODEL: ORDER PARAMETER AND SPECTRAL DENSITY
The Lagrangian density of the Thirring model in threedimensional Euclidean space-time has the form and contains a vector-vector interaction built from N flavors ψ 1 , . . . , ψ N . In the present work ψ a (or ψ) always denotes a 4-component reducible spinor. The hermitean matrices Γ µ with µ = 1, 2, 3 form a 4-dimensional reducible representation of the Clifford algebra. After introducing a Hubbard-Stratonovich auxiliary vector field v µ , a subsequent integration over the fermion fields leads to the partition function (see [1] for more details) Here we used that the determinant for N flavors is just the N 'th power of the determinant for 1-flavor with Dirac operator We introduced a chirality-breaking fermion mass which is needed to control our lattice Monte-Carlo simulations in the chirally broken phase. The eigenvalues of iD come in pairs (λ+i m, −λ+i m) such that the fermion determinant is real and positive, This means that the effective action in (2) is real or that the (massive or massless) Thirring model with N reducible flavors has no sign problem. Hence in the well-known auxiliary field formulation the model can be simulated by Monte-Carlo methods on a space-time lattice. At this point we observe that N is just a parameter that can be varied continuously. In the present work we will focus on N 1 and thus consider lattice models which continuously extrapolate to N = 1.0 from below. The so defined models have no parity anomaly for any real N . The massless Thirring model with N reducible flavors is invariant under the discrete Z 2 parity transformation as well as global U(2N ) chiral transformations. These symmetries, together with the discrete C and T symmetries, are well explained in [24]. A technical problem here is that on a finite lattice the condensates vanish in the massless case exactly for every vector field configuration and a careful extrapolation to vanishing fermion mass is difficult.
For performance reasons, we simulate the theory in a 2-component irreducible representation of the Clifford algebra. A convenient reducible representation is and the corresponding Dirac operator (3) reads At this point we change the fermionic variables, ψ a → Γ 45 ψ a ,ψ a →ψ a , a = 1, . . . , N, such that the Dirac operator / D (acting on two-component irreducible spinors) enters D with the same sign 2 , i.e. that D in (6) is replaced by The effective action in (2) takes the form As order parameter for chiral symmetry we use the chiral condensate where the insertion Σ 45 originates from the change of variables in (7). Using translational invariance it can be written as We see here that only the Dirac operator / D of one irreducible flavor -introduced in (6) -enters the expression for the partition function and chiral condensate of N reducible flavors. Note that the condensate defined in (10) is real and positive, Σ = |Σ|. In terms of the spectral density ρ v of the irreducible Dirac-operator in a fixed auxiliary field, defined by the condensate (11) can be written as where the non-negative expectation valueρ(E) is calculated with the effective action, The last relation follows from charge conjugation symmetry which implies ρ v = ρ −v and is explained in appendix A. In the limit m → 0 equation (13) gives rise to a variant of the celebrated Banks-Casher relation [32]. It relates the low end of the spectral density of the irreducible operator i / D to the chiral condensate of the reducible models. In passing we note that -because of parity-symmetry -the would-be order parameter of parity ∝ ψ a ψ a is identically zero for all reducible models. This means that there is no spontaneous breaking of parity. Finally, we must emphasize that varying the number of reducible flavors N continuously between 1 2 and 1 as described above is not equivalent to varying the number of irreducible flavors N ir between 1 and 2. There are several reasons for this difference: First and rather technically, the N ir = 1 model suffers from a severe sign problem and can only be simulated in an interesting dual formulation [1], in contrast to the reducible model with N = 1 2 , which has no sign problem. Second and more important, for N ir = 1 the Z 2 parity symmetry can be broken (by the anomaly and/or spontaneously) while parity is never broken for the reducible systems.
In the Thirring models with N ir = 2 and N = 1 the global U(2) chiral symmetry can be broken to U(1)×U (1) in which case we should see two massless Goldstone bosons in the particle spectrum. Finally we note, that the interpolating models with N / ∈ N/2 probably do not describe local quantum field theories. But this problem will not invalidate the reasoning in the present work.

III. CHIRAL CONDENSATE
We performed simulations with chiral SLAC-fermions on lattices L × (L − 1) 2 in the range L = 6 . . . 24. To control and stabilize our simulations, we chose a mass proportional to the inverse lattice size, with small dimensionless parameter m 0 . Note that for any fixed value of m 0 one recovers the massless Thirring model in the infinite volume limit L → ∞. For small λ the condensate vanishes due to the (annoying but well-known) large lattice artifacts in the strong coupling regime, [1,27].
In order to determine the critical flavor number, we investigate the maximum Σ max of the λ-dependent chiral condensate Σ for different flavour numbers N and lattice sizes L. The maximum of the condensate is well-motivated since it clearly signals the breaking of chiral symmetry. The obtained results fully comply with those obtained with the alternative method based on the susceptibility of the interaction term in a later section. Figure 2 shows the dependence of Σ max on the mass parameter m 0 for three different lattice volumes and for N = 0.70. For a fixed m 0 (with Compton wave-length much smaller than the lattice size) the chiral condensate increases with increasing lattice volume. Performing the infinite volume limit -which includes the m → 0 limit for every m 0 > 0 -we conclude that for N = 0.7 chiral symmetry is spontaneously broken. Actually, in most of our simulations we choose m 0 = 0.1, which is a good compromise between good chiral properties, simulation performance and small finite volume effects. The results for the maximal condensate Σ max as function of N (for m 0 = 0.1) is depicted in Figure 3. For a fixed lattice volume, the condensate increases with decreasing flavor number. For a fixed N 0.8 the maximal condensate increases with increasing lattice volume and we conclude that chiral symmetry is broken for these N . We compared with the results obtained with m 0 = 0.04 and obtained a comparable outcome. But for this smaller mass finite size effects are less suppressed. The region above N = 0.8 is magnified in Figure 4. Above N = 0.95 the condensate decreases with increasing volume and one concludes that chiral symmetry remains unbroken in this regime. Unfortunately, the lattices are not sufficiently large to permit a reliable extrapolation to infinite volume for all values of N under consideration. That was only achieved for the flavor numbers below 0.75 and above 1.00. Three examples are depicted in Figure 5. Since we introduced a mass, we expect a finite size scaling law of the form [33] Σ max (L) = ae −bL + Σ max (∞) (16)   for which the optimal fit-parameters in the fits depicted in Figure 5 are listed in Table I. 3 In the broken phase with small N (e.g. 0.75) this extrapolation works well. Also for N = 1.0 the exponential function (16) fits the data well and points to a vanishing condensate in the infinite volume limit. For values of N in between the data becomes basically flat due to large finite size effects -in some cases they are even non-monotonic -which renders an extrapolation unreliable. However, for every finite volume we find that the maximal chiral condensate exhibits a turning point around N ≈ 0.8 where the chiral condensate is bending upwards, see Figure 3. This bending is caused by finite size effects and the explicit breaking of chiral symmetry by the fermion mass term. The data points to the left of this turning point are well-described by the scaling law with parameters a, N c , β given in Table II. In particular we can read off the critical flavour number and conclude, that there is no spontaneous symmetry breaking above In the following sections we will substantiate the result (18) with other methods. Note that our lattice volumes are not large enough to extract a reliable value for the critical exponent β. But since our main focus is on the critical flavor number, which does not suffer from finite size effects, we did not further increase the lattice volume to obtain a  [30] more accurate value for β. The critical exponent β has been calculated previously with the functional renormalization group (FRG), with Dyson-Schwinger equations (DSE) and with Monte-Carlo simulation with staggered fermions (MC). We compiled some results with references in Table III. We see that the predictions for the critical exponent β depend much on the non-perturbative method in use. The quoted values cannot be easily compared among themselves and with our results in Table II. For example, with staggered fermions one may simulate another universality class. We intend to find a better value of β with chiral fermions on larger lattices in the future.
For the smaller mass parameter m 0 = 0.04 we obtain qualitatively the same data. However, the ill-conditioned fermion determinant forbids a more detailed study for this (and smaller) masses.

IV. SPECTRAL DENSITY
As explained above, the chiral properties of the theory can be extracted from the spectral density ρ v (E) of the massless irreducible Dirac operator introduced in (12) and the average spectral densityρ(E) defined in (14). Ifρ(E) in the neighborhood of E = 0 remains small with increasing volume, then chiral symmetry is realized. On the contrary, if it increases, then chiral symmetry is broken. Figure 6 shows the spectral density for N = 0.80 on different lattice sizes. Close to the origin, the density clearly builds up with increasing lattice volume and one concludes that chiral symmetry is broken. For the larger flavor number N = 1.00 we observe the opposite behavior, see Figure 7: Close to the origin, the density remains small for all lattice sizes. Again we conclude that for N = 1.00 chiral symmetry is unbroken.  (14) for N = 1.00 (symmetric phase) for different lattice sizes. The shaded regions indicate the uncertainties.

V. GOLDSTONE SPECTRUM
Next we investigate the meson spectrum of the N -flavor theory. There are two scalar and two pseudoscalar mesons with vanishing angular momentum and the corresponding interpolating operators are O Γ = 1 N aψ a Γψ a , where Γ is the identity matrix or one of the three matrices iΓ 4 , iΓ 5 and Γ 45 in (5). Since all reducible flavors contribute equally to O Γ , we may set N = 1 in these bilinears. Thus we choose the operator basis which are the zero-momentum projections of and where ψ represents one of the N reducible flavors. For example, σ 1 ⊗ σ 0 swaps the two irreducible spinors which make up the reducible 4-component spinor. Note that the expectation value of O 3 (x) is twice the chiral condensate. In our simulations, we measure the correlation matrix with elements where ∆ is the propagator for 4-component fermions in a fixed auxiliary field v µ , The expectation values in (21) are calculated with S eff and traces are taken in spinor and flavor space. By exploiting parity invariance we prove in appendix A that the correlation matrix is diagonal. It is most conveniently expressed in terms of the parity odd and parity even terms in the decomposition where The diagonal elements of (C ab ) -these are the eigenvalues -read C 0 (t) = 4 x ,y where x = (t, x ) and y = (0, y ). If chiral symmetry is spontaneously broken according to C 1 = C 11 and C 2 = C 22 should describe massless particles.
If chiral symmetry is not broken, we expect that the four (pseudo)scalars arrange in a singlet and a triplet of SU(2) ⊂ U(2). In particular, the state belonging to the interpolating operators O 1 , O 2 and O 3 should form a triplet. In the corresponding subspace the correlation matrix has eigenvalues C 1 , C 2 and C 3 . Indeed, in the symmetric phase we have B xy = 0 for m → 0 and these 3 eigenvalues become identical, In Figure 8 we show the (pseudo)scalar spectrum in the symmetric phase at N = 1.00 for two different lattice volumes 11×24 and 15×24 and a residual mass m = 0.004. The correlation functions C 1 , C 2 and C 3 for both spatial volumes lie almost on top of each other -the splitting originates from the explicit breaking by the mass term -while C 0 decays faster. The lines represent fits with a sum of two cosh-functions for the ground and excited state. The fitted masses are given in Table IV. For both the ground and excited multiplet we find three almost identical masses and a larger one. Within statistical uncertainties and taking into account finite volume effects, the results are compatible with two multiplets of massive mesons in the symmetric phase.
In the broken phase at N = 0.80, see Figure 9, the correlation functions C 1 = C 2 and C 3 differ significantly compared to the correlators in the symmetric phase. While the masses of C 1 = C 2 are almost volume independent, the ground state mass of C 3 shows stronger finite volume effects. The correlation function C 0 is compatible with zero for all temporal extents t which is explained by a corresponding correlation length not much bigger than the  In an interesting recent work by Simon Hands with bulk domain wall fermions (DWF) on a 12 3 -lattice (and L s up to 40) the meson correlators of the N = 1 model have been calculated as well [3]. Whereas an earlier simulation with surface DWF on a 12 2 × 24 lattice (but L s only up to 16) showed no sign of a chiral phase transition for N = 2 [34], the new results for N = 1 with DWF signal a Goldstone spectrum expected from a U(2)→ U(1) ⊗ U(1) breaking. This means that for 1 flavor the prediction of DWF are still in conflict with our findings. based on different approaches that there exists a welldefined continuum limit, corresponding to a UV-stable fixed point of the renormalization group (RG) [15,16,20,23,35]. To find the continuum theory at the transition to the artifact phase at strong bare couplings -see [1] for details -seems unlikely since this transition only exists in a discretized setup. However, already in the quoted work we have spotted signals of another transition in the intermediate coupling regime. In this section we will argue, that such a transition indeed exists for N ≥ N c and probably is continuous. In our earlier work we did not further analyze this feature, mainly since scanning the phase diagram of a fermion theory on lattices of different sizes is rather expensive. For the same reason we do not aim at a detailed finite size analysis in the present work. But we do simulations on lattices with different sizes to see the qualitative behavior of the susceptibility related to the four-Fermi term in the Lagrangian. Actually, the similarly accurate number for N c is extracted by spotting the merging of the newly discovered transition with the lattice artifact transition 4 .
As tracer for the transition we will consider the second derivative of the partition function with respect to the coupling λ. As discussed in detail in our previous paper [1], the partition function's first derivative can be associated with the lattice filling factor k as follows, which is (up to an additive constant) proportional to the expectation value of the four-Fermi term (ψγ µ ψ) 2 . Roughly speaking, k is the average fraction of lattice sites at which an interaction takes place. From this interpretation the following properties (established in [1]) are comprehensible: The filling factor vanishes in the weak coupling regime (large λ), it monotonically approaches 1 when approaching the lattice artifact phase at strong coupling and its derivative exhibits a dip (inverted peak) at this transition. All of these features are clearly seen in Figure 10 and Figure 11. What has not been discussed before is the fact that for not too small N the derivative ∂ λ k, which is related to the second derivative of the thermodynamic potential ln Z and hence the susceptibility of the operator (ψΓ µ ψ) 2 , shows a second dip for intermediate values of λ. This was already visible at the edge of Figure 4 in [1]. Since the direct computation of ∂ λ k as 8-point function would be rather expensive computationally, we instead use the numerical derivative of k to calculate the susceptibility. But conventional finite-difference approximations of the λ-derivative will greatly amplify the noise present in our data. There are many methods to regularize the differentiation process (regression, smoothing, filtering, variation denoising). In our analysis we used a variation denoising method (and checked it with the conventional approach). More details can be found in appendix B.
Examples of ∂ λ k at three different flavour numbers N are depicted in Figure 11. One recognizes two qualitatively different behaviors. For small N (N = 0.7 in Figure 11) ∂ λ k has one distinct minimum which -as discussed before -signals the transition into the lattice artifact phase. For sufficiently large N (N = 1.0 in Figure 11) two distinct minima are clearly visible. A comparison with the data presented in Figure 10 reveals that the minimum to the left (strong coupling) signals the transition into the lattice artifact phase. The second peak at intermediate coupling has not been discussed before and the following discussion makes clear, that it belongs to an interaction This value matches the putative critical flavor number N c in (18) pretty well. Our explanation of this only seemingly surprising equality is the following: the ubiquitous lattice artifact phase at strong coupling does not describe any properties of the continuum Thirring model. Only in the physical phase at weaker couplings can we hope to construct a continuum theory when approaching a critical point or critical line of second order transitions. For sufficiently small N the perfect candidate for this transition is the chirality breaking transition discussed previously. Indeed, the line where the condensate is maximal is always to the right of both minima of ∂ λ k, see singularity. The second order chiral transition turns into a second order transition without order parameter. The order of a phase transition is related to the dependence of the peak susceptibility on the size of the system. For N = 1.0 this behavior of k and ∂ λ k is depicted in Fig. 13 while the depth of the minima and the position of the interaction-driven transition are shown in Fig. 14 and  15 respectively. In both figures, one can see that finite size effects are significant for L < 14. Above that, at the artifact transition the susceptibility ∝ ∂ λ k is almost independent on the volume as expected for a first order transition. On the other hand, at the interaction-driven transition the susceptibility increases with increasing volume roughly according to [36] ln χ max (L) = a ln L + b as expected for a second order transition. Fitting this behavior to the maximal susceptibilities of both transitions (for L ≥ 14), we find The theory of finite-size scaling furthermore predicts that the coupling λ max (L), where χ(L) peaks, approaches the critical coupling in the thermodynamic limit λ c as [36] λ max (L) = λ c (1 − cL −1/ν ).
We observe that at the interaction-driven transition this scaling law reproduces the date well on the whole range of lattices sizes in use with the values λ c = 0.526(2), c = 1(2) · 10 2 , ν = 0.26 (6). (34) The data and corresponding fit are depicted in Figure  15. Since the linear coefficient a from above should obey a = γ/ν − d (with dimension d = 3), we find the following rough estimate for the critical exponent γ, The aim of our rather crude finite size analysis is not to find accurate critical exponents but rather to study the order of the interaction-driven transition. Our results clearly suggest that it is a second order transition at which the correlation length should diverge. Most likely this continuous transition is not associated with any symmetry breaking, since the term (ψΓ µ ψ) 2 is already part of the Thirring-Lagrangian (1). Such transitions without change of symmetry are common in condensed matter physics and they are sometimes called iso-symmetric. In solid state physics such transitions are structural and are related to discontinuous changes of the cell volume and cell parameters and thus indicate a first-order transition. But continuous transitions without symmetry breaking are also possible in which case we prefer the name interaction-driven transitions. For example, a continuous transition without symmetry-breaking bilinear fermion condensate -triggered by a four-fermion interaction term -has been reported previously in SU(4)invariant four-fermi models in three dimensions. These models are similar to the Thirring model considered in the present work. Numerical simulations with staggered fermions, the fermion bag method, hybrid Monte Carlo and quantum Monte Carlo revealed actually an interesting phase structure [37][38][39][40]: while increasing the strength of the four-Fermi term in these models. They exhibit a continuous quantum phase transition from a weakly coupled massless phase (a gapless Dirac semimetal) to a massive (fully gapped Mott insulator) phase without condensing any fermion bilinear operator. It could very well be that a similar mechanism is at work in the Thirring models, although a bilinear condensate is not forbidden by symmetry arguments as it is in the SU(4)-invariant models.
Although the transition is probably not associated with a change of symmetry there could be an order parameter and the filling factor k is a possible candidate. From Figure 13 one might conjecture that for weak coupling (to the right of the interaction-driven transition) k approaches 0 in the infinite volume limit. This would imply that only the phase between the peaks describes an interacting four-Fermi theory. Actually, we can prove that k = 1 in the strong coupling expansion, see [1], but so far we could not show that k = 0 in a weak coupling expansion.
To summarize: we conjecture that the critical number N c separating systems with and without chiral symmetry breaking and the number N t where the two phase transition lines come close or meet should be identified.
As function of λ the chiral condensate decreases rapidly to the left of the peak -towards the strong coupling region -and decreases slowly towards the weak coupling region. But how can we explain that on a finite lattice the chiral condensate is maximal for N ≥ N t just to the right of the interaction-driven transition line in Fig. 12 and only vanishes in the infinite volume limit? Without SSB there are two sources for a non-zero condensate: the explicit symmetry breaking by the fermion mass and fluctuations. After (11) we have argued that any v µ -configuration adds a non-negative number to the condensate. Near a second order transition the fluctuations are large and on a finite system these large fluctuations drive the condensate away from zero. This explains, why the condensate peaks near the interaction-driven transition line. On the other hand, near a first order transition to the lattice artifact phase the fluctuations do not necessarily grow and we do not expect a fluctuation driven condensate. This is the reason, why for N ≥ N t (in which case the first-order and second-order lines are well-separated) the condensate is small near the artifact transition line and does not depend much on the volume. Then we would expect, that the condensate just to the right of the artifact line is a better approximation to the condensate in the thermodynamic limit. In Fig. 16 we plotted for every N the maximum of the chiral condensate and its value in the proximity -actually just to the right -of the lattice artifact transition line 5 . We see that for N < N t the chiral condensate Σ prox follows the old fit in Fig. 3 (with the form (17) and the parameters from Table II). This is expected since N t ≈ N c . Nevertheless, it further substantiates our claim that the condensate Σ prox is a better approximation to the chiral condensate at infinite volume as compared to the maximal condensate since fluctuations, which drive the condensate away from the infinite volume result, are suppressed.

VII. CONCLUSIONS
In the present work we have re-analyzed the longstanding problem about the critical flavor number in the three-dimensional (reducible) Thirring models. We used chiral SLAC fermions to have full control over the chiral properties of the model. In this formulation the chiral and parity symmetry are manifest and no fine tuning is required. We reformulated the model such that the number of reducible flavors N becomes a continuous parameter -offering the possibility of determining precisely when spontaneous symmetry breaking ceases to exist. We calculated the maximum of the chiral condensate, the spectral density and the spectrum of scalar and pseudo-scalar particles as function of the flavor number N 5 Just to the right means three ticks (in the fixed λ-grid) to the right. For comparison we applied this rule to the maxima of the condensate and the artificial phase transition line. The points where Σprox in Fig. 16 have been measured are depicted in Fig. 12.
between 0.5 and 1.0. As a result we find a critical flavor number N c = 0.80 (4) .
In particular, we spotted two Goldstone bosons only for N ≤ N c . Since a non-integer value of N probably does not describe a local quantum field theory (and in particular no Thirring model), we conclude that there is no spontaneous symmetry breaking in all reducible Thirring models.
With an elaborate and expensive scan of the susceptibility related to the interaction term (ψΓ µ ψ) 2 as function of the coupling λ and the number of flavors N we spotted -besides the (probably first order) ubiquitous lattice artifact transition -a (probably second order) transition for all N greater than We gave several arguments why N c and N t should be identified. We expect that for an arbitrary number of flavors N there exists a continuum limit: for every N ≤ N c = N t there exists a QFT with spontaneous breaking of chiral symmetry and for every N ≥ N c there exists a chirally invariant QFT. This fully agrees with our previous findings in [1] and other recent lattice Monte-Carlo approaches: The domain wall formulation in [3] supports our claim that there is no spontaneous symmetry breaking for N > 1. Since there are still major technical issues to be studied in the domain wall formulation -such as the discrepancy between the bulk and surface formulation and the additional L s → ∞ extrapolation -the conclusion for the N = 1 case is only preliminary. But it seems to disagree with the results in the present analysis with SLAC fermions and in [1].
In parallel to the present work L. Dabelow, H. Gies and B. Knorr investigated reducible Gross-Neveu-Thirring models in three dimensions with FRG methods by admitting momentum dependent vertices in the flow equation for the scale dependent effective action [31]. Their new estimate for N c (obtained with their most strict criterion) is compatible with ours.
We would like to stress that our results are not in contradiction with those in [1,41], where a breaking of parity symmetry in models with an odd number of irreducible flavors has been reported. The irreducible models are very different from the parity invariant reducible models studied in the present work and in other more recent publications on the three-dimensional Thirring model.
Besides the question about the precise value of N c we witness a convergence of recent results obtained with sophisticated functional methods and lattice simulations based on chiral fermions. So the question arises why earlier attempts with staggered fermions failed to predict an acceptable value for N c ? It has already been pointed out in [1,3], and we would like to stress it once more, that the failure of staggered fermions to find the correct symmetry (or even universality class) and phase structure of 3-dimensional four-Fermi theories away from weak coupling, is probably also responsible for the mismatch between DMF and staggered fermion results near a conformal fixed point in 3 + 1 dimensional non-Abelian gauge theory [42]. For strongly coupled (fermion) systems we should be careful to implement all global internal symmetries in any discretization.
Simulations of fermion systems are rather time consuming and an elaborate finite size analysis could not be accomplished in the present work. For example, to decide about the order of the chiral phase transition below N c and the interaction-driven phase transition above N c requires further extensive studies. Even more demanding would it be to extract critical exponents of interest to decide about the universality class of the system at criticality. This would allow for a comparison with recent results obtained wih functional methods. We hope to report on further progress in these directions in the near future.
where tr F denotes the trace in flavor space, and in addition Recall, that the correlation matrix (21) involved suitable traces over the spinor indices as well, tr = tr D tr F . Next we study the transformation of the Greenfunction under parity. Since the eigenmodes change according to (A2) and the eigenvalues swap signs, we have ∆ ± (x, y, v) = −γ 3 ∆ ∓ (x,ỹ,ṽ)γ 3 , In the last step we used, that the two γ 3 in the conjugation (A6) chancel under the trace over Dirac indices and that summing over all x is the same as summing over allx .
Averaging with the parity-invariant effective action over the auxiliary field results into Similarly one obtains x ,y tr D A xx B yy = 0 = x ,y It follows that the correlation matrix C(t) is diagonal with eigenvalues C a (t) given in (25). Finally note that i tr D B xx is just the chiral condensate Σ.

Appendix B: Numerical Differentiation
While numerical differentiation of smooth data is easily done by discrete derivative stencils, non-smooth and particularly noisy data is hard to differentiate numerically. This is seen in Figure 11 where the markers show the result of applying the stencil to the rather smooth looking data of Figure 10. Particularly, in the interesting regime around N = 0.8 such a numerical differentiation is basically useless because of the large noise. Another approach, which we will use in the following, is total-variation (TV) regularized differentiation [43]. It reformulates the problem as a global optimization problem such that the minimum of the functional F (u) = I(u) − (k − k(λ 0 )) + αR(u) is assumed for an approximation u ≈ ∂ λ k. Here I(u) is an (appropriate discrete) integration operation and · an appropriate norm such that k is obtained from integrating its derivative u. Afterwards a regulator term R can be added to smooth the minimizing solution u. 6 While one can clearly see the smoothing behavior of this approach in Figure 11, the important information about the peak (location) is not distorted compared to the naive scheme. We always cross-checked that the TV result was plausible within the naive scheme; however, we cannot assign a pointwise uncertainty to the TV result due to the global procedure for obtaining it.  .
The coefficients α i and β i depend on the degree n of the approximation, on the power of the inverse fermion matrix k and on the spectral range of M . Details on the RHMC algorithm can be found in [44]. In this way we are able to perform lattice simulations for any rational flavour number. To speed up the simulations, we use different approximations in the HMC trajectory and the metropolis acceptance step. For most of our simulations we use p = 4 pseudo-fermions and a degree of the approximation n HMC , n acc = 10, 25. The inverse of the shifted fermion matrix in the rational approximation is computed by a multi-mass conjugate gradient (CG) solver. During the CG iterations, we have to apply the SLAC operator to a pseudo-fermion field. Here, we make use of a special property of the SLAC derivative: It is diagonal in momentum space and we obtain where the sum is over all lattice momenta p. Instead of using a three-dimensional (parallelized) Fourier transformation, we apply one-dimensional Fourier transformations that are computed in parallel. Although there is communication overhead, this method is on small lattices far more efficient than a three-dimensional Fourier transformation.