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

TL;DR

This study investigates the Thirring model in 2+1 dimensions, finding no evidence of chiral symmetry breaking for any integer number of fermion flavors, and identifies a critical flavor number around 0.8.

Contribution

The paper provides an accurate determination of the critical flavor number in the Thirring model using analytic continuation, challenging previous claims of symmetry breaking.

Findings

Critical flavor number N^c ≈ 0.80(4).

No chiral symmetry breaking for N ≥ 1.

Evidence of a new phase transition for N > N^c.

Abstract

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.

Figures

Tables

Table 1. Table 1: Fit parameter for the infinite volume extrapolation of the maximal chiral condensate for three values of N 𝑁 N .

$N$	$a$	$b$	${\Sigma }_{\max }\left(\infty \right)$
$0.75$	$-5.9\left(1\right)$	$0.10\left(2\right)$	$4.21\left(4\right)$
$0.9$	$-1.0\left(3\right)$	$0.20\left(9\right)$	$1.29\left(6\right)$
$1.0$	$1\left(1\right)$	$0.03\left(4\right)$	$0\left(1\right)$

Table 2. Table 2: Fit parameters for fitting the critical behavior.

$L$	$N$	$a$	$\beta$	$N^{\text{c}}\left(L\right)$
$8$	$[0.50,0.75]$	$0.245\left(4\right)$	$0.6\left(1\right)$	$0.82\left(4\right)$
$12$	$[0.50,0.80]$	$0.266\left(3\right)$	$0.31\left(2\right)$	$0.815\left(3\right)$
$16$	$[0.50,0.79]$	$0.276\left(4\right)$	$0.23\left(2\right)$	$0.801\left(4\right)$

Table 3. Table 3: Compilation of various results for the critical exponent β 𝛽 \beta from the literature. The abbreviations are explained in the main text.

Method	$\beta$	Ref.
FRG	0.44	Janssen and Gies (2012)
DSE	1	Kondo (1995)
MC (staggered)	0.37	Christofi et al. (2007)

Table 4. Table 4: Meson masses in the symmetric phase with N = 1.00 𝑁 1.00 N=1.00 and the broken phase with N = 0.80 𝑁 0.80 N=0.80 for two different spatial lattices L s = 11 subscript 𝐿 s 11 L_{\text{s}}=11 and L s = 15 subscript 𝐿 s 15 L_{\text{s}}=15 .

$C$	$m\left(11\right)$	$m\left(15\right)$	$m^{\ast }\left(11\right)$	$m^{\ast }\left(15\right)$	$N$
$C_0$	$0.21\left(2\right)$	$0.21\left(2\right)$	$1.27\left(6\right)$	$1.22\left(7\right)$	$1.00$
$C_{1,2}$	$0.134\left(3\right)$	$0.128\left(2\right)$	$1.03\left(5\right)$	$1.02\left(3\right)$	$1.00$
$C_3$	$0.138\left(2\right)$	$0.131\left(2\right)$	$1.08\left(4\right)$	$0.98\left(3\right)$	$1.00$
$C_0$	$-$	$-$	$-$	$-$	$0.80$
$C_{1,2}$	$0.103\left(2\right)$	$0.095\left(3\right)$	$1.04\left(12\right)$	$0.93\left(17\right)$	$0.80$
$C_3$	$0.109\left(4\right)$	$0.127\left(7\right)$	$0.81\left(7\right)$	$0.81\left(10\right)$	$0.80$

Equations

$$\mathcal{L}=\sum _{a=1}^N{\bar{\psi }}_a\mathrm{i}{\Gamma }^{\mu }{\partial }_{\mu }{\psi }_a-\frac{g^2}{2N}{j}^{\mu }{j}_{\mu },j^{\mu }=\sum _{a=1}^N{\bar{\psi }}_a{\Gamma }^{\mu }{\psi }_a,$$

$$Z$$

$$S_{\mathrm{eff}}$$

$$D={\Gamma }^{\mu }\left({\partial }_{\mu }-\mathrm{i}{v}_{\mu }\right)+m.$$

$$\det (\mathrm{i}D)={\left(\det D^{†}D\right)}^{\frac{1}{2}}\ge 0.$$

$${\Gamma }^{\mu }=$$

$${\Gamma }_5=$$

$$D=\left(\begin{array}{cc}/D+m& 0\\ 0& -/D+m\end{array}\right),/D={\gamma }^{\mu }({\partial }_{\mu }-\mathrm{i}{v}_{\mu }).$$

$${\psi }_a\to {\Gamma }_{45}{\psi }_a,{\bar{\psi }}_a\to {\bar{\psi }}_a,a=1,\dots ,N,$$

$$D=\left(\begin{array}{cc}/D+m& 0\\ 0& /D-m\end{array}\right),/D={\gamma }^{\mu }({\partial }_{\mu }-\mathrm{i}{v}_{\mu }).$$

$$S_{\mathrm{eff}}(v)=\lambda\int\mathrm{d}^{3}x\,v^{\mu}v_{\mu}-N\,\ln\det\big{(}m^{2}-\not{D}^{2}\big{)}\,.$$

$$\Sigma =\frac{\mathrm{i}}{2N}\sum _a⟨{\bar{\psi }}_a{\Sigma }_{45}{\psi }_a⟩,$$

$$\Sigma =\frac{1}{V}\frac{1}{Z}\int \mathcal{D}{v}^{\mu }\mathrm{tr}\left(\frac{m}{m^2-/D^2}\right)e^{-S_{\mathrm{eff}}(v)}.$$

$$\mathrm{tr}f(\mathrm{i}/D)={\int }_{-\infty }^{\infty }\mathrm{d}Ef(E){\rho }_v(E),$$

$$\Sigma =\frac{2m}{V}{\int }_0^{\infty }\frac{\mathrm{d}E}{E^2+m^2}\bar{\rho }(E),$$

$$\bar{\rho }(E)=\frac{1}{Z}\int \mathcal{D}{v}^{\mu }{e}^{-S_{\mathrm{eff}}(v)}{\rho }_v(E)=\bar{\rho }(-E).$$

$$m=\frac{m_0}{L},$$

$${\Sigma }_{\max }(L)=a{e}^{-bL}+{\Sigma }_{\max }(\infty )$$

$${\Sigma }_{\max }(N)=a(N^{\text{c}}-N{)}^{\beta }$$

$$N^{\text{c}}=0.80(4).$$

$${\mathcal{O}}_a(t)=\sum _{\text{boldmathe}x}{\mathcal{O}}_a(t,\text{boldmathe}x),a=0,1,2,3,$$

$${\mathcal{O}}_a(x)=\bar{\psi }(x)({\sigma }_a\otimes {\sigma }_0)\psi (x)={\mathcal{O}}_a(t,\text{boldmathe}x),$$

$$C_{ab}(t)$$

$$\displaystyle=\sum_{\boldmathe{x},\boldmathe{y}}\big{\langle}\operatorname{tr}{\sigma_{a}\Delta_{xx}}\operatorname{tr}{\sigma_{b}\Delta_{yy}}-\operatorname{tr}(\sigma_{a}\Delta_{xy}\sigma_{b}\Delta_{yx})\big{\rangle}$$

$$-\sum _{\text{boldmathe}x,\text{boldmathe}y}⟨\mathrm{tr}{\sigma }_a{\Delta }_{xx}⟩⟨\mathrm{tr}{\sigma }_b{\Delta }_{yy}⟩,$$

$$\Delta =\frac{1}{D}=\left(\begin{array}{cc}{\Delta }_{+}& 0\\ 0& {\Delta }_{-}\end{array}\right),{\Delta }_{\pm }=⟨x|\frac{1}{\mathrm{i}/D\pm \mathrm{i}m}|y⟩.$$

$$\Delta ={\sigma }_0\otimes A+{\sigma }_3\otimes B,$$

$$A=\frac{\mathrm{i}/D}{m^2-/D^2},\mathrm{i}B=\frac{m}{m^2-/D^2}.$$

$$\displaystyle\begin{split}C_{0}(t)&=4\sum_{\boldmathe{x},\boldmathe{y}}\big{\langle}\operatorname{tr}_{\mathrm{\hskip 0.85358ptD}}A_{xx}\operatorname{tr}_{\mathrm{\hskip 0.85358ptD}}A_{yy}\big{\rangle}\\ &-2\sum_{\boldmathe{x},\boldmathe{y}}\big{\langle}\operatorname{tr}_{\mathrm{\hskip 0.85358ptD}}(A_{xy}A_{yx}+B_{xy}B_{yx})\big{\rangle}\,,\\ C_{1}(t)&=\;C_{2}(t)=2\sum_{\boldmathe{x},\boldmathe{y}}\big{\langle}\operatorname{tr}_{\mathrm{\hskip 0.85358ptD}}(B_{xy}B_{yx}-A_{xy}A_{yx})\big{\rangle}\,,\\ C_{3}(t)&=4\sum_{\boldmathe{x},\boldmathe{y}}\big{(}\big{\langle}\operatorname{tr}_{\mathrm{\hskip 0.85358ptD}}B_{xx}\operatorname{tr}_{\mathrm{\hskip 0.85358ptD}}B_{yy}\big{\rangle}+\Sigma^{2}\big{)}\\ &-2\sum_{\boldmathe{x},\boldmathe{y}}\big{\langle}\operatorname{tr}_{\mathrm{\hskip 0.85358ptD}}(A_{xy}A_{yx}+B_{xy}B_{yx})\big{\rangle}\,,\end{split}$$

$$\text{U}(2)\to \text{U}(1)\otimes \text{U}(1)$$

$$C_0(t)$$

$$C_1(t)$$

$$k=-\frac{\lambda}{NV}\frac{\partial\ln Z(\lambda)}{\partial\lambda}+\mathrm{c}=\frac{1}{4N\lambda}\,\big{\langle}j^{\mu}j_{\mu}\big{\rangle}+\mathrm{c}\,.$$

Full text

I Introduction

The Thirring model Thirring (1958) in $2$ space-time dimensions is integrable and in the massless limit even soluble Klaiber (1968); Sachs and Wipf (1996). The model in $3$ space-time dimensions is of interest for various reasons, e.g. its close relationship to QED3 Hands (1995); Itoh et al. (1995); Gusynin and Pyatkovskiy (2016); Kotikov and Teber (2016) or its relevance in solid state physics, where it describes low-energy electronic properties of materials like graphene Semenoff (1984); Hands and Strouthos (2008) or high-temperature superconductors Herbut (2002); Franz et al. (2002). In $3$ dimensions the model is perturbatively non-renormalizable but can be renormalized in the limit of large flavor numbers $N$ Parisi (1975); Gawedzki and Kupiainen (1985); Rosenstein et al. (1991); Hands (1995). Thus it provides a simple realization Braun et al. (2011) of the concept of asymptotic safety Weinberg (1976). In the large-$N$ limit one finds an unbroken U$(2N)$ symmetry for every coupling strength. On the other hand, in the limit $N=\nicefrac{{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 couplings111More precisely, the Thirring model with $1$ irreducible $2$-component Fermion flavor is the same as the Gross-Neveu model with $1$ irreducible Fermion flavor.. We conclude that the Thirring model exhibits no chiral phase transition for large $N$ but shows a second order phase transition at $N=\nicefrac{{1}}{{2}}$. The question about the critical flavor number $N^{\text{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^{\text{c}}=2$ to $N^{\text{c}}=\infty$ Gomes et al. (1991); Itoh et al. (1995); Sugiura (1997); Kondo (1995); Gies and Janssen (2010); Janssen and Gies (2012); Kim and Kim (1996); Del Debbio and Hands (1996); Del Debbio et al. (1997); Del Debbio and Hands (1999); Hands and Lucini (1999); Christofi et al. (2007), more recent lattice studies with chiral fermions favor smaller values of $N^{\text{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\gtrapprox 1$ Wellegehausen et al. (2017). 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 Hands (2019). The results support our finding that $N^{\text{c}}$ is smaller than hitherto believed with the notable difference that he interprets his data as an evidence for $1<N^{\text{c}}<2$. In a recent explorative functional renormalization group (FRG) study with momentum-dependent couplings and the pseudo-spectral method the critical behavior of four-fermion theories Dabelow et al. (2019) has been reconsidered. While a precise estimate for $N^{\text{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\geq 1.0$ but also verify the proposed SSB $\mathrm{U}(2)\to\mathrm{U}(1)\otimes\mathrm{U}(1)$. We conclude that indeed there is a critical flavor number $N^{\text{c}}\approx 0.80$ which is considerably smaller than $1.0$. A similarly accurate value for $N^{\text{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 interaction driven phase transition for all models with flavor numbers $N\geq N^{\mathrm{t}}=0.78(4)$. We argue that $N^{\mathrm{t}}$ should be identified with $N^{\text{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,\dots$ 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 Wellegehausen et al. (2017), 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 $\propto(\bar{\psi}\Gamma_{\mu}\psi)^{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 three-dimensional Euclidean space-time has the form

[TABLE]

and contains a vector-vector interaction built from $N$ flavors $\psi_{1},\dots,\psi_{N}$. In the present work $\psi_{a}$ (or $\psi$) always denotes a $4$-component reducible spinor. The hermitean matrices $\Gamma_{\mu}$ with $\mu=1,2,3$ form a $4$-dimensional reducible representation of the Clifford algebra. After introducing a Hubbard-Stratonovich auxiliary vector field $v^{\mu}$, a subsequent integration over the fermion fields leads to the partition function (see Wellegehausen et al. (2017) for more details)

[TABLE]

Here we used that the determinant for $N$ flavors is just the $N$’th power of the determinant for $1$-flavor with Dirac operator

[TABLE]

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 $\mathrm{i}D$ come in pairs $(\lambda+\mathrm{i}\,m,-\lambda+\mathrm{i}\,m)$ such that the fermion determinant is real and positive,

[TABLE]

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\lessapprox 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 $\mathbbm{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 Janssen and Gies (2012). 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

[TABLE]

and the corresponding Dirac operator (3) reads

[TABLE]

At this point we change the fermionic variables,

[TABLE]

such that the Dirac operator $\not{D}$ (acting on two-component irreducible spinors) enters $D$ with the same sign222If $\not{D}$ would have a different sign for the two irreducible flavors, then Euclidean correlators could violate positivity constraints., i.e. that $D$ in (6) is replaced by

[TABLE]

The effective action in (2) takes the form

[TABLE]

As order parameter for chiral symmetry we use the chiral condensate

[TABLE]

where the insertion $\Sigma_{45}$ originates from the change of variables in (7). Using translational invariance it can be written as

[TABLE]

We see here that only the Dirac operator $\not{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 is real and positive, $\Sigma=|\Sigma|$. In terms of the spectral density $\rho_{v}$ of the irreducible Dirac-operator in a fixed auxiliary field, defined by

[TABLE]

the condensate (11) can be written as

[TABLE]

where the non-negative expectation value $\bar{\rho}(E)$ is calculated with the effective action,

[TABLE]

The last relation follows from charge conjugation symmetry which implies $\rho_{v}=\rho_{-v}$ and is explained in appendix A. In the limit $m\to 0$ equation (13) gives rise to a variant of the celebrated Banks-Casher relation Banks and Casher (1980). It relates the low end of the spectral density of the irreducible operator $\mathrm{i}\not{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 $\propto\langle\bar{\psi}_{a}\psi_{a}\rangle$ 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 $\frac{1}{2}$ and $1$ as described above is not equivalent to varying the number of irreducible flavors $N_{\text{ir}}$ between $1$ and $2$. There are several reasons for this difference: First and rather technically, the $N_{\text{ir}}=1$ model suffers from a severe sign problem and can only be simulated in an interesting dual formulation Wellegehausen et al. (2017), in contrast to the reducible model with $N=\frac{1}{2}$, which has no sign problem. Second and more important, for $N_{\text{ir}}=1$ the ${\mathbb{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_{\text{ir}}=2$ and $N=1$ the global U$(2)$ chiral symmetry can be broken to U$(1)\times$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\notin{\mathbb{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\times(L-1)^{2}$ in the range $L=6\dots 24$. To control and stabilize our simulations, we chose a mass proportional to the inverse lattice size,

[TABLE]

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\to\infty$.

Figure 1 shows the surface plot of the chiral condensate for $\lambda=0.25\dots 0.60$ and $N=0.5\dots 1.1$ on a $16\times 15^{2}$ lattice with $m_{0}=0.1$.

For $N=0.5$ there exists a broken phase with non-zero chiral condensate $\Sigma$. With increasing flavor number the chiral condensate falls off very quickly. For small $\lambda$ the condensate vanishes due to the (annoying but well-known) large lattice artifacts in the strong coupling regime, Del Debbio et al. (1997); Wellegehausen et al. (2017).

In order to determine the critical flavor number, we investigate the maximum $\Sigma_{\mathrm{max}}$ of the $\lambda$-dependent chiral condensate $\Sigma$ for different flavor 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 $\Sigma_{\mathrm{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\to 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 $\Sigma_{\mathrm{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\lessapprox 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 Neuberger (1989)

[TABLE]

for which the optimal fit-parameters in the fits depicted in Figure 5 are listed in Table 1.333In all the following fits, $L=14$ was excluded from fitting because there were some incurable problems with the thermalization. Also, $\Sigma(\infty)$ was constrained to be positive, since the exact condensate has this property.

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\approx 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

[TABLE]

with parameters $a,N^{\text{c}},\beta$ given in Table 2. In particular we can read off the critical flavor number and conclude, that there is no spontaneous symmetry breaking above

[TABLE]

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 $\beta$. 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 more accurate value for $\beta$.

The critical exponent $\beta$ 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 3.

We see that the predictions for the critical exponent $\beta$ depend much on the non-perturbative method in use. But the quoted values cannot be easily compared among themselves and with our results in Table 2. For example, with staggered fermions one may simulate another universality class. We intend to find a better value of $\beta$ 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 $\rho_{v}(E)$ of the massless irreducible Dirac operator introduced in (12) and the average spectral density $\bar{\rho}(E)$ defined in (14). If $\bar{\rho}(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.

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 $\mathcal{O}_{\Gamma}=\frac{1}{N}\sum_{a}\bar{\psi}_{a}\Gamma\psi_{a}$, where $\Gamma$ is the identity matrix or one of the three matrices $\mathrm{i}\Gamma_{4},\mathrm{i}\Gamma_{5}$ and $\Gamma_{45}$ in (5). Since all reducible flavors contribute equally to $\mathcal{O}_{\Gamma}$, we may set $N=1$ in these bilinears. Thus we choose the operator basis

[TABLE]

which are the zero-momentum projections of

[TABLE]

and where $\psi$ represents one of the $N$ reducible flavors. For example, $\sigma_{1}\otimes\sigma_{0}$ swaps the two irreducible spinors which make up the reducible $4$-component spinor. Note that the expectation value of $\mathcal{O}_{3}(x)$ is twice the chiral condensate. In our simulations, we measure the correlation matrix with elements

[TABLE]

where $\Delta$ is the propagator for $4$-component fermions in a fixed auxiliary field $v_{\mu}$,

[TABLE]

The expectation values in (21) are calculated with $S_{\mathrm{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

[TABLE]

where

[TABLE]

The diagonal elements of $(C_{ab})$ – these are the eigenvalues – read

[TABLE]

where $x=(t,\boldmathe{x})$ and $y=(0,\boldmathe{y})$.

If chiral symmetry is spontaneously broken according to

[TABLE]

two of the four mesons should become massless Goldstone bosons while the other two remain massive. More precisely, since the interpolating operators $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$ correspond to the Goldstone bosons, their correlators $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)\subset$ U$(2)$. In particular, the state belonging to the interpolating operators $\mathcal{O}_{1},\mathcal{O}_{2}$ and $\mathcal{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\to 0$ and these $3$ eigenvalues become identical,

[TABLE]

In Figure 8 we show the (pseudo)scalar spectrum in the symmetric phase at $N=1.00$ for two different lattice volumes $11\times 24$ and $15\times 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 4. 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 lattice constant. Thus the simulation results are compatible with the existence of two massless Goldstone bosons and two massive excited states with different masses. We conclude, that the chiral U$(2)$ symmetry is indeed broken to U$(1)\otimes\text{U}(1)$ for $N=0.80$.

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 Hands (2019). Whereas an earlier simulation with surface DWF on a $12^{2}\times 24$ lattice (but $L_{s}$ only up to $16$) showed no sign of a chiral phase transition for $N=2$ Hands (2016), the new results for $N=1$ with DWF signal a Goldstone spectrum expected from a U(2)$\to$ U(1)$\,\otimes$ U(1) breaking. This means that for $1$ flavor the prediction of DWF are in conflict with our findings.

VI Susceptibility

If for $N\geq N^{\text{c}}$ the lattice artifact phase transition were the only phase transition then one could hardly imagine how to construct an interacting QFT in the continuum limit. And there are convincing arguments based on different approaches that there exists a well-defined continuum limit, corresponding to a UV-stable fixed point of the renormalization group (RG) Parisi (1975); Hikami and Muta (1977); Gomes et al. (1991); Gawedzki and Kupiainen (1985); Gies and Janssen (2010). To find the continuum theory at the transition to the artifact phase at strong bare couplings – see Wellegehausen et al. (2017) 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\geq N^{\text{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^{\text{c}}$ is extracted by spotting the merging of the newly discovered transition with the lattice artifact transition 444We cannot be sure that the two transitions meet. But we will see that they come close..

As tracer for the transition we will consider the second derivative of the partition function with respect to the coupling $\lambda$. As discussed in detail in our previous paper Wellegehausen et al. (2017), the partition function’s first derivative can be associated with the lattice filling factor $k$ as follows,

[TABLE]

and thus is (up to a $\lambda$-independent additive constant c) proportional to the expectation value of the four-Fermi interaction term $(\bar{\psi}\gamma^{\mu}\psi)^{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 Wellegehausen et al. (2017)) are comprehensible: The filling factor vanishes in the weak coupling regime (large $\lambda$), it monotonically approaches $1$ when approaching the lattice artifact phase at strong coupling and its derivative

[TABLE]

exhibits a dip at this transition, since the susceptibility of $j^{\mu}j_{\mu}$ – this is just the sum over the first argument of the connected two-point function on the right hand side – peaks at the 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 $\partial_{\lambda}k$ shows a second dip corresponding to a second peak of the susceptibility of the $4$-Fermi term at intermediate values of $\lambda$. This was already visible at the edge of Figure 4 in Wellegehausen et al. (2017). Since the direct computation of $\partial_{\lambda}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 $\lambda$-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 compared it with the conventional approach). More details can be found in appendix B.

Examples of $\partial_{\lambda}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) $\partial_{\lambda}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 (inverted) peak at intermediate coupling has not been discussed before and the following discussion makes clear that it belongs to an interaction driven transition at intermediate couplings. For flavor numbers near $N=0.8$ the two minima merge and the artifact phase transition line comes close to or meets the line of interaction driven-phase transitions.

The numerically extracted peak positions of the susceptibility for $N$ between $0.7$ and $1.0$ are shown in Figure 12. At first glance one can see that for $N<0.76$ there is only one phase transition and for $N>0.82$ there are two transitions. More accurately, on a lattice with $L=16$ the two transition lines get close or meet at the triple point at

[TABLE]

This value matches the putative critical flavor number $N^{\text{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 $\partial_{\lambda}k$, see Fig. 12. Hence for $N\geq N^{\mathrm{t}}$ the maximum is to the right of the interaction driven phase transition line and we have seen that this maximal condensate decreases with increasing lattice size. The only plausible explanation is the following: for small $N\leq N^{\text{c}}$ there is, beside the artifact transition at strong coupling a second order chiral phase transition at intermediate coupling. The solid line in Fig. 12 indicates the position of the maximal chiral susceptibility (obtained as numerical derivative $\chi_{\Sigma}=-\partial_{\lambda}\Sigma$) where the chiral phase transition on sufficiently large lattices happens. At $\lambda\approx 0.42$ and $N=N^{\text{c}}$ the chiral phase transition line comes close to the artifact transition line and ceases to exist since in the thermodynamic limit there is no non-vanishing condensate for $N\geq N^{\text{c}}$. Instead a (probably) second order phase transition line emerges for $N\geq N^{\mathrm{t}}\approx N^{\text{c}}$ where the derivative $\partial_{\lambda}k$ develops a singularity. The second order chiral transition with order parameter turns into a second order interaction-driven 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 $\partial_{\lambda}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 $\propto\partial_{\lambda}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 Janke (2008)

[TABLE]

as expected for a second order transition, where $\chi\propto-\partial_{\lambda}k+\mathrm{const}$ is the susceptibility of the $4$-Fermi term, see (29). Such a linear dependence is actually seen in the double-logarithmic plot in Figure 14 for lattices with $L\geq 14$. However, the data points for lattices with $L\geq 14$ are too noisy and we cannot extract reliable values for the fit parameters $a$ and $b$. The theory of finite-size scaling also predicts that the coupling $\lambda_{\mathrm{max}}(L)$, where $\chi(L)$ peaks, approaches the critical coupling in the thermodynamic limit $\lambda_{\mathrm{c}}$ as Janke (2008)

[TABLE]

We observe that at the interaction-driven transition this scaling law reproduces the data reasonably well. The extracted value

[TABLE]

is rather stable – it does not change much if the fits are for lattice with $L\geq 8,\,L\geq 10$ or $L\geq 12$, see Figure 15. But the extracted values for $\nu$ vary considerably and cannot be trusted. The lattices considered are just not big enough to determine the critical exponent $\nu$. But the aim of our very crude finite size analysis is not to calculate critical exponents but rather to study the order of the interaction-driven transition. Our preliminary results suggest that it is a second order transition with infinite correlation length.

Most likely this continuous transition is not associated with any symmetry breaking, since the term $(\bar{\psi}\Gamma_{\mu}\psi)^{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 Ayyar and Chandrasekharan (2015); Catterall (2016); He et al. (2016); Ayyar and Chandrasekharan (2016): The systems 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 exist an order parameter. 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 [math] 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 Wellegehausen et al. (2017), but so far we could not show that $k=0$ in a weak coupling expansion.

To summarize: we conjecture that the critical number $N^{\text{c}}$ separating systems with and without chiral symmetry breaking and the number $N^{\mathrm{t}}$ where the two phase transition lines come close or meet should be identified.

But how can we explain that on a finite lattice the chiral condensate is maximal for $N\geq N^{\mathrm{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^{\mu}$-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 and the chiral susceptibility555For $N\geq N^{\text{c}}$ there is still a peak of the susceptibility shown as dashed line in Fig. 12. $\chi_{\Sigma}$ peak 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\geq N^{\mathrm{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 line666Just to the right means three ticks (in the fixed $\lambda$-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 $\Sigma_{\mathrm{prox}}$ in Fig. 16 have been measured are depicted in Fig. 12.. We see that for $N<N^{\mathrm{t}}$ the chiral condensate $\Sigma_{\mathrm{prox}}$ follows the old fit in Fig. 3 (with the form (17) and the parameters from Table 2). This is expected since $N^{\mathrm{t}}\approx N^{\text{c}}$. Nevertheless, it further substantiates our claim that the condensate $\Sigma_{\mathrm{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 long-standing 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$ between $0.5$ and $1.0$. As a result we find a critical flavor number

[TABLE]

In particular, we spotted two Goldstone bosons only for $N\leq N^{\text{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 $(\bar{\psi}\Gamma_{\mu}\psi)^{2}$ as function of the coupling $\lambda$ 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

[TABLE]

We gave several arguments why $N^{\text{c}}$ and $N^{\mathrm{t}}$ should be identified. Thus we expect that for an arbitrary number of flavors there exists a continuous phase transition: for every $N\leq N^{\text{c}}=N^{\mathrm{t}}$ there is a transition with spontaneous breaking of chiral symmetry and for every $N\geq N^{\text{c}}$ there exists a transition without spontaneous breaking of chiral symmetry. But since $N^{\text{c}}<1$ only the latter transition can be used to construct continuum Thirring models with $N=1,2,3,\dots$.

The result (34) improves the result in Wellegehausen et al. (2017) and results of lattice Monte-Carlo simulations with domain wall fermions in Hands (2019). The latter support 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}\to\infty$ 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 Wellegehausen et al. (2017).

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 Dabelow et al. (2019). Their new estimate for $N^{\text{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 Schmidt et al. (2016); Wellegehausen et al. (2017), 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^{\text{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^{\text{c}}$? It has already been pointed out in Wellegehausen et al. (2017); Hands (2019), 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 Hasenfratz et al. (2017). 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 really decide about the order of the interaction-driven phase transition above $N^{\text{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 with functional methods. We hope to report on further progress in these directions in the near future.

Acknowledgments

We are grateful to Shailesh Chandrasekharan, Lukas Janssen, Simon Hands, Rajamani Narayanan and in particular Holger Gies for helpful discussions and comments. We thank Holger Gies for communicating to us their FRG-results prior to publication. Daniel Schmidt contributed considerably at an earlier stage of this project. Our simulations were performed on the HPC cluster at the University of Jena.

Appendix A Parity

We choose the parity operation $x\to\tilde{x}=(x_{1},x_{2},L-x_{3})$, where $L$ is the extend of the box in $3$-direction. The auxiliary vector field transforms as

[TABLE]

and a $2$-component spinor field $\chi$ as

[TABLE]

Now it follows at once that if $\chi$ is an eigenfunction of $\mathrm{i}\not{D}_{v}$ with eigenvalue $\lambda$, then $\tilde{\chi}$ is an eigenfunction of $\mathrm{i}\not{D}_{\tilde{v}}$ with eigenvalue $-\lambda$.

A.1 Spectral density

The invariance under parity implies particular properties of the spectral density defined in (12) and the fermion Greenfunction $\Delta$ of the reducible Diracoperator in (8). For example, we conclude

[TABLE]

Since the space-time integrals over $v_{\mu}^{2}$ and $\tilde{v}_{\mu}^{2}$ are equal and in addition only the square of $\not{D}$ enters the effective action, we see that the latter is parity invariant,

[TABLE]

Since the averaging over the auxiliary field is done with the parity invariant factor $e^{-S_{\mathrm{eff}}}$ we end up with the relation (14) which states, that the averaged spectral density $\bar{\rho}(E)$ is an even function of the spectral parameter.

A.2 Fermion Greenfunction

The fermion Greenfunction of the reducible system is

[TABLE]

where $\not{D}$ belongs to the irreducible system. The Greenfunction is a linear combination of $\sigma_{0}$ and $\sigma_{3}$, see (23) and (24). It follows that

[TABLE]

where $\operatorname{tr}_{\mathrm{\hskip 0.85358ptF}}$ denotes the trace in flavor space, and in addition

[TABLE]

Recall, that the correlation matrix (21) involved suitable traces over the spinor indices as well, $\operatorname{tr}=\operatorname{tr}_{\mathrm{\hskip 0.85358ptD}}\operatorname{tr}_{\mathrm{\hskip 0.85358ptF}}$. Next we study the transformation of the Greenfunction under parity. Since the eigenmodes change according to (37) and the eigenvalues swap signs, we have

[TABLE]

which in turn implies

[TABLE]

It follows, for example, that

[TABLE]

In the last step we used, that the two $\gamma_{3}$ in the conjugation (41) chancel under the trace over Dirac indices and that summing over all $\boldmathe{x}$ is the same as summing over all $\tilde{\boldmathe{x}}$. Averaging with the parity-invariant effective action over the auxiliary field results into

[TABLE]

Similarly one obtains

[TABLE]

It follows that the correlation matrix $C(t)$ is diagonal with eigenvalues $C_{a}(t)$ given in (25). Finally note that $\,\mathrm{i}\operatorname{tr}_{\mathrm{\hskip 0.85358ptD}}\langle B_{xx}\rangle$ is just the chiral condensate $\Sigma$.

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

[TABLE]

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 Chartrand (2011). It reformulates the problem as a global optimization problem such that the minimum of the functional

[TABLE]

is assumed for an approximation $u\approx\partial_{\lambda}k$. Here $I(u)$ is an (appropriate discrete) integration operation and $\|\cdot\|$ 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$.777Details of the implementation can be found in Chartrand (2011). We used the python translation from https://github.com/stur86/tvregdiff of their MATLAB code with minor modifications. 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.

Appendix C Monte-Carlo simulations

For performance reasons, the Monte-Carlo simulations have been performed in the two-component irreducible representation. In the HMC algorithm we compute the fermion determinant by introducing $p$ pseudo-fermions $\phi,\phi^{\dagger}$

[TABLE]

where the function $r(A,n)$ is a rational approximation of the inverse fermion matrix

[TABLE]

and

[TABLE]

The coefficients $\alpha_{i}$ and $\beta_{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 Clark and Kennedy (2004). In this way we are able to perform lattice simulations for any rational flavor 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_{\text{HMC}},n_{\text{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

[TABLE]

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.