On the Critical Flavor Number in the 2+1$d$ Thirring Model

The Thirring model in 2+1 spacetime dimensions, in which $N$ flavors of relativistic fermion interact via a contact interaction between conserved fermion currents, is studied using lattice field theory simulations employing domain wall fermions, which furnish the correct U(2N) global symmetry in the limit that the wall separation $L_s\to\infty$. Attention is focussed on the issue of spontaneous symmetry breakdown via a non-vanishing fermion bilinear condensate $\langle\bar\psi\psi\rangle\not=0$. Results from quenched simulations are presented demonstrating that a non-zero condensate does indeed form over a range of couplings, provided simulation results are first extrapolated to the $L_s\to\infty$ limit. Next, results from simulations with $N=1$ using an RHMC algorithm demonstrate that U(2) symmetry is unbroken at weak coupling but plausibly broken at strong coupling. Correlators of mesons with spin zero are consistent with the Goldstone spectrum expected from U(2)$\to$U(1)$\otimes$U(1). We infer the existence of a symmetry-breaking phase transition at some finite coupling, and combine this with previous simulation results to deduce that the critical number of flavors for the existence of a quantum critical point in the Thirring model satisfies $01$.


Introduction
While the study of theories of relativistic fermions moving in the plane, ie. in 2+1 spacetime dimensions, has received a fillip over the past decade as a result of developments in the condensed matter physics of layered systems, the underlying quantum field theories continue to be of considerable theoretical interest in their own right. This paper concerns the d = 2 + 1 Thirring model, describing N flavors of interacting fermion, whose Lagrangian density reads Here µ = 0, 1, 2, and a key feature is that the fields ψ,ψ lie in reducible spinor representations, so that Dirac matrices γ µ are 4 × 4. This enables the definition of a mass term mψψ which is parity-invariant, where it is convenient to define a discrete parity transformation in terms of inversion of all three spacetime axes.
Once m = 0 only (2) remain as symmetries; ie. dynamical fermion mass generation corresponds to a breaking pattern U(2N) →U(N)⊗U(N). The question of whether bilinear condensation takes place is inherently non-perturbative, and has been first studied using truncated Schwinger-Dyson equations [2,4,5]. While details are somewhat scheme-dependent, the picture emerging is that symmetry breaking is possible for sufficiently large g 2 and sufficiently small N, and indeed there exists a critical flavor number N c such that no symmetry breaking occurs at any coupling for N > N c . More recently the Functional Renormalisation Group (FRG) has also been applied [6]. The main focus of this work has been to establish the existence of UV-stable RG fixed points g 2 * (N), such that an interacting field theory exists at all scales in the limit g 2 → g 2 * . On the assumption that the symmetry-breaking transition is second-order, it seems reasonable to identify this quantum critical point (QCP) with the critical g 2 c (N), which exists for N < N c . The identification of N c is thus an important ingredient in the search for novel QCPs.
It is natural to apply lattice field theory methods to the problem, and indeed this has been tried by several groups over the years [7] - [12]. A common feature of all these approaches is the use of staggered lattice fermions in 2+1d. The conclusion of [11], employing simulation studies in the effective strong-coupling limit, is that N c = 6.6(1), and that the critical exponent δ characterising the response of the order parameter to an explicit symmetry-breaking mass at criticality has value δ(N c ) ≈ 7. Away from the strong-coupling limit the value of δ is found to be rather sensitive to N. This is to be compared with the Schwinger-Dyson predictions N c ≃ 4.32, δ(N c ) = 1 [4]. In summary, the staggered Thirring model exhibits a non-trivial phase diagram in the (N, g 2 ) plane with some interesting features.
At the level of the lattice action, massless staggered fermions in 2+1d have a manifest U(N stag )⊗U(N stag ) global symmetry, distinct from the U(2N) of eqn. (1), which is broken to U(N stag ) by a fermion mass. In a weak-coupling continuum limit, N stag staggered fermions are known to describe N = 2N stag continuum flavors, with an eventual recovery of U(2N) at long wavelengths [13]. Near a QCP, however, the story may be different. Indeed, a study of the N = 2 staggered model using a fermion bag algorithm [12], which permits simulations directly in the massless limit, found critical exponents compatible with those of the Gross-Neveu model [14], an unxepected result since on the face of it the two models have distinct Lagrangians, different symmetries, and completely different 1/N expansions, with symmetry breaking due to bilinear condensation predicted in the large-N limit in the GN case, and expected to persist for all N. Instead, the results of [12,14] imply the two models lie in the same RG basin of attraction. Indeed, when written purely in terms of four-point interactions between staggered fermion fields spread over the vertices of elementary cubes, the only difference between Thirring and GN is an extra body-diagonal coupling in the latter case [14].
The mismatch between theoretical expectation and results from the staggered Thirring model has motivated us to consider alternative lattice fermion formulations with the potential to capture the resquisite symmetries more faithfully. Specifically, we have developed both analytical and numerical insight into how U(2N) symmetry is manifested in 2+1d using Domain Wall Fermions (DWF) [15,16], which will be reviewed in the next section. Next, in Ref. [17] we applied DWF in exploratory simulations of both GN and Thirring models with N = 2. The main conclusions of that work were that in the GN model there appears to be no obstruction to studying symmetry breaking via bilinear condensation and the resulting QCP; good qualitative agreement was found with the analytical expectations of the large-N approach. However, no evidence was found for symmetry breaking in the Thirring model with fixed L s = 16, implying N c < 2 in contradiction to the staggered fermion model. Meanwhile, the Jena group has applied another U(2N)-invariant formulation, the SLAC fermion, and found no symmetry breaking in the Thirring model all the way down to N = 1 [18]. We also note in this context the recent reported mismatch between DWF and staggered fermion results near a conformal fixed point in 3+1d non-abelian gauge theory, which potentially springs from the same failure of the staggered action to capture the correct symmetries away from weak coupling [19]. This paper extends the study of spontaneous symmetry breaking in the Thirring model with DWF to N = 1. The formulation is reviewed below in Sec. 2. Of the two possible ways to introduce the Thirring interaction discussed in [17], we will focus almost exclusively on the variant in which the auxiliary vector field is located uniformly throughout the bulk, ie. stressing its resemblence to an abelian gauge potential. Simulating just a single DWF flavor requires the RHMC algorithm, as argued in Appendix A; details of the implementation are also given in Sec. 2, along with results of an initial survey at fixed L s = 8, permitting a comparison with the N = 2 data from [17]. Since to date U(2N) symmetry breaking has not been observed with DWF, as a warm-up Sec. 3 presents results from a study of the quenched Thirring model, to see what broken symmetry might look like. We will see that extrapolation to the L s → ∞ limit is crucial at the stronger couplings examined, and introduce an exponential Ansatz which empirically works well. There does indeed appear to be a range of couplings where U(2N) symmetry is broken. Results from an exploratory study of a quenched model with the auxiliary formulated just on the domain walls are also presented. Finally in Sec. 4 we present a detailed study of the N = 1 model on 12 3 at four representative g 2 over a range of masses m, allowing L s to vary between 8 and 40 (and in some cases 48) to facilitate for the first time a controlled L s → ∞ extrapolation. Comparison data taken for N = 2, and for N = 1 on 16 3 are also presented. In addition to the bilinear condensate ψ ψ , results for correlators of spin-0 mesons will be given, including the channel with quantum numbers of the would-be Goldstone bosons. We will argue that at the strongest coupling examined, close to the effective strong coupling limit, the most plausible explanation of the data is that U(2) symmetry is spontaneously broken. By contrast, the L s -extrapolated data for the N = 2 model in the effective strong coupling limit is consistent with unbroken symmetry, implying 1 < N c < 2. We summarise and outline plans for future work in Sec. 5.

Formulation and Implementation
The most straightforward way to simulate the Thirring model with orthodox techniques is via introduction of a vector auxiliary field A µ (x). The continuum Lagrangian density is then written In this form the similarity to an abelian gauge theory is manifest, and it is clear the Thirring model inherits the same U(2N) global symmetry (2,3). The bosonic action violates the gauge symmetry, however, though this can be remedied by introduction of a Stückelberg scalar leading to a "hidden local symmetry" [4]. Eqn. (4) includes a mass term mψψ, which for reducible spinor representations in 2+1d can be shown to be invariant under parity inversion. However, it is not unique; as outlined in [15] there are three possible parity-invariant mass terms: m hψ ψ; im 3ψ γ 3 ψ; im 5ψ γ 5 ψ; in Euclidean metric the first is hermitian while the two "twisted" forms are antihermitian. Due to their equivalence under U(2N) rotations, and the absence of chiral anomalies in 2+1d, they are physically indistinguishable. In this work fermions are studied on a 2+1+1d lattice using the DWF formulation. The action is written [17] where Ψ(x, s),Ψ(x, s) are defined on the 2+1+1d lattice. S bos is the action for the auxiliary boson fields A µ (x) defined on links µ = 0, 1, 2, and is an obvious generalisation of the gaussian term in (4). For convenience throughout we will use lattice units with a ≡ 1, but note here that the dimensionless combination is g −2 a. The coordinates x, y denote sites in 2+1d, and s running along the third direction x 3 takes values 1, . . . , L s . The fermion kinetic action is defined D W (M) x,y is the 2+1d Wilson operator with M the domain wall height: while D 3 governs hopping along x 3 : The factors (1−δ s ′ ,1/Ls ) implement open boundary conditions at domain walls located at s = 1, L s , while the projectors P ± ≡ 1 2 (1 ± γ 3 ) also appear in the identification of target fermions ψ(x),ψ(x) defined as 2+1d fields localised on the domain walls at s = 1, L s : The relations (10) are a major ingredient in the physical interpretation of DWF, but for now we stress they should be regarded as assumptions. They permit a definition of the mass term m a S a in (7), using (5) with a = h, 3, 5. It is easily checked that m h , m 3 couple fields defined on opposite walls, while m 5 couples fields on the same wall. In Ref. [16] it was shown that in the limit L s → ∞ the DWF operator M is equivalent to an overlap operator constructed to satisfy 2+1d generalisations of the Ginsparg-Wilson relations, and [15] demonstrated recovery of U(2N) symmetry in a weakly-coupled theory, quenched non-compact QED 3 , in the same limit. Specifically, for finite L s the three condensates related by U(2N) satisfy: all residuals decaying exponentially with L s with hierarchy δ h ≫ ǫ h ≫ ǫ 3 , ǫ 5 . The dominant residual δ h is defined by the imaginary component of the 3-condensate: and is thus measurable even when, as in this work, the action m 3 S 3 is used.
To complete the specification of the Thirring model with DWF we need the fermionauxiliary interaction S int . Two variants were introduced in Ref. [17]. The surface formulation has the link fields A µ linearly coupled to point-split fermion bilinears defined on the walls: A similar approach has been adopted in DWF studies of another 2+1d theory of interacting fermions, the Gross-Neveu model [20,17]. Because the interaction is defined only at the walls, (15) brings the technical advantage that the Pauli-Villars determinant needed to formally recover the correct fermion measure as L s → ∞ does not depend on A µ , and hence need not be simulated, making calculations with (15) relatively inexpensive. The bulk formulation emphasises the resemblence of the vector auxiliary to a gauge field, defining a linear interaction between the vector bilinear current and an s-independent A µ throughout the bulk: (16) Operationally, the fermion operator with S int = S bulk resembles that of a gauge theory with connection U µ = (1 + iA µ ); in other words, the link field is no longer constrained to be unitary. This choice is not unique -other lattice approaches to the Thirring model use unitary link fields [8,21] -but for N > 1 it ensures there are no fermion interactions higher than four-point once the auxiliary is integrated out [7]. However, in the same work it was shown that this regularisation fails to preserve the transversity of the vacuum polarisation correction to the A-propagator, leading to an additive renormalisation of g −2 in the large-N expansion. The consequent uncertainty in identifying the strongcoupling limit has been explored using staggered fermions in [11], where the pragmatic approach of identifying the physical strong coupling limit g 2 R → ∞ with the point where ψ ψ(g 2 ) has a maximum was found to yield a plausible equation of state. This point will be further discussed below Fig. 1.
Most of the results presented in the paper are obtained using the bulk formulation, and it will be shown in Sec. 4.3 why this is the preferred option. Appendix A derives some relations for the fermion determinant in the bulk approach, motivating the use of the RHMC algorithm [22] for simulation with N = 1. For even N, the HMC algorithm outlined in [17] is sufficient. The pseudofermion action used in RHMC is where subscripts on M † M denote non-vanishing mass terms. We choose the domain wall mass M = 1. The components with m h = 1 describe Pauli-Villars fields needed to cancel bulk contributions to the determinant and ensure coincidence with the correct overlap operator as L s → ∞ [16]. The fractional matrix powers needed to compute (17) are estimated by means of a rational approximation The required coefficents α, β are calculated with the Remez algorithm using the implementation available at [23]. They were chosen such that over a spectral range (0.0001, 50) (which accommodates the upper limit (2(2 + 1 + 1) − M) 2 obtained in the free-field limit of (6)), |r p (x) − x p | is less than 10 −6 for matrices needed during guidance and 10 −13 for those needed in the Hamiltonian calculations required for the acceptance step of the algorithm. This appears to be a conservative requirement for the systems studied to date, and translates into N pf = 12 (guidance) and N pf = 25 (acceptance). The need for further refinement cannot be ruled out for future studies of critical systems. The partial fraction expansion (18) is efficiently calculated using the multi-shift procedure described in [24], which in turn requires the use of a hermitian Lanczos solver such as described in [25]. For the systems we have examined, particularly as L s is made large, maintaining orthonormality of the Lanczos vectors generated at each successive iteration requires double precision arithmetic; on the same systems the conjugate gradient algorithm used in measurement routines runs happily in single precision. For an evaluation of x = (M † M) p Φ, the convergence criterion adopted is where ρ i [25] is a real variable parametrising the magnitude of the latest increment to the solution vector x i (where x = α 0 Φ + α i x i ), and ε = 10 −6 (guidance) and 10 −9 (acceptance). Finally, it should be noted that the matrix inversions are numerically demanding, especially as the coupling becomes strong, possibly as a consequence of the non-unitarity of the link felds U µ . For instance, on the largest 16 3 × 40 volume studied, at the strongest coupling g −2 = 0.3 and the smallest mass m = 0.01, the Lanczos solver in the Hamiltonian calculation requires roughly 11000 iterations to achieve convergence. The conjugate gradient solver operating on stochastic noise sources in the measurement routine on the same system requires roughly 4700 iterations.   For orientation, in the remainder of this section we present the main features of the model based on simulations with fixed finite L s . All results are taken using the m 3 S 3 mass term, and henceforth for convenience when there is no possibility of confusion we will often write the associated condensate ψ iγ 3 ψ as ψ ψ . Fig. 1 shows ψ ψ vs g −2 for m = 0.01 on 12 3 × L s , for both surface and bulk formulations with N = 1 (L s = 8) and N = 2 (L s = 16). The N = 2 results obtained using the HMC algorithm were first presented in [17]. In all cases the condensate increases as the coupling is increased from weak to strong, until it reaches a maximum in the region g −2 ≈ 0.2 -0.3. This non-monotonic behaviour maximum is also observed in simulations using both staggered [7,11] and SLAC [26] fermions, and is associated with strong-coupling artifacts possibly due to the non-transversity of the auxiliary propagator discussed above; following [11] we will identify the maximum with the approximate location of the effective strong coupling limit, and focus our attention on the weak-coupling side of this maximum. With the vertical scale chosen to accommodate the N = 1 bulk data, the condensates obtained with the surface formulation in this region are very small and show little dependence on coupling. The bulk formulation yields larger condensates, but the most striking feature of Fig. 1 is the sharp rise in the bulk condensate for g −2 < 0.6 for N = 1; it is already apparent that the tendency for fermions and antifermions to pair is more significant here than for any case previously examined. Fig. 2 shows the auxiliary action density N 2V g −2 xµ A 2 µ (x) vs. g −2 on the same systems. This measurable offers an interesting diagnostic of the UV properties of the different model approaches. First note that in the continuum, a comparison of ∂ ln Z/∂g 2 obtained using the original action (1) and with the bosonised form (4), following a change in functional integration variables, results in the following identity for the boson action: Assuming smooth behaviour of the expectation of the square of the fermion current on the RHS of (20), we therefore expect departures from the free-field value 3 2 to increase monotonically with g 2 , and Fig. 2 shows this is indeed the case for the surface model with N = 1, 2. For the bulk model the corresponding relation contains terms of the formΨγ µ Ψ(s)Ψγ µ Ψ(s ′ ),Ψγ µ Ψ(s)Φ † γ µ Φ(s ′ ), with s, s ′ = 1, . . . , L s and Φ, Φ † are Pauli-Villars fields, ie. there are contributions from bulk fields whose interpretation is not as transparent. In fact, Fig. 2 shows the correction to the free-field result has the opposite sign except at the very strongest coupling.
The contrast between surface and bulk formulations was already noted in [17]. Here we note that N = 1, 2 yield very similar results for the surface formulation, but for N = 1 there is a marked contrast in the bulk results once g −2 < ∼ 0.5, again hinting at interesting strong coupling behaviour. Fig. 2 also plots N = 1 bulk data from 12 3 × 40 showing small but significant disparities; this is a reminder of the importance of seeking the L s → ∞ limit of all observables in the DWF approach, in particular the bilinear condensate. Sec. 3 presents a first investigation in this direction in the quenched limit N = 0, and enables us to address the question what does spontaneous symmetry breaking due to bilinear condensation look like with DWF? .

Results in the Quenched Limit N = 0
The quenched theory with N = 0 is technically very simple to explore; one simply performs fermionic measurements using the operator M on field configurations inexpensively generated using the gaussian auxiliary action. Unlike gauge theory, there is no theoretical expectation that the results have any relevance to the full theory; this is best understood via the auxiliary propagator S A (x) in the large-N expansion, which at strong coupling decays as |x| −2 as a result of vacuum polarisation corrections [3], but which remains a contact S A (x) ∼ δ d (x) in the quenched limit. Fig. 3 compares condensate data obtained with the bulk formulation for N = 0, 1, 2 for m = 0.01 with L s = 16; for N > 0 the spacetime volume is 12 3 , but the low computational cost enabled the quenched study on 16 3 . Compared to Fig. 1 the vertical scale has been extended to accommodate the quenched data: the hierarchy ψ ψ(N = 0) ≫ ψ ψ(N = 1) ≫ ψ ψ(N = 2) is as expected; the low eigenvalues of the effective Dirac operator responsible for the condensate signal via the Banks-Casher relation also suppress the determinant in the path integral measure.  To explore the L s → ∞ limit we performed a systematic study of the bilinear condensate ψ ψ(g 2 , m) on 16 3 × L s with L s = 8, . . . , 40, g −2 = 0.2, 0.3, . . . , 1.0 and m = 0.01, 0.02, . . . , 0.05. Each boson configuration, separated by 100 HMC trajectories, was analysed using 10 stochastic noise vectors located on either wall. 25000 trajectories were studied for L s = 8, 16, and 5000 for L s = 24, 32, 40. Results for ψ ψ(L s ) with m = 0.05 and varying g −2 are shown in Fig. 4, and for varying m at g −2 = 0.4, 0.8 in Fig. 5. It is evident that finite-L s corrections are significant, and increase in importance as the coupling grows. We have modelled them using the notation of (12) as follows: The resulting three parameter fits are plotted as dashed lines in Figs. 4,5. The exponential form (21) works well across the dataset, but the asymptotic value ψ ψ ∞ becomes poorly constrained as m → 0 resulting in large uncertainties in this limit. Also note that the strongest coupling g −2 = 0.2 looks to be an outlier in both Figs. 3 and 4, reflecting the probable influence of strong coupling artifacts.  Results for ψ ψ ∞ obtained using fits to (21) are plotted for various g −2 as a function of mass m in Fig. 6. The curves are for the most part remarkably m-independent. The m → 0 limit of the strong coupling data at g −2 = 0.2 and possibly 0.3 are affected by artifacts as discussed above, and the large errorbars in the same limit for g −2 > ∼ 0.8 reflect poorly constrained fit parameters associated with the lack of curvature seen in Fig. 5. It is difficult to say anything definitive in either case. However. Fig. 6 supports a coupling window g −2 ∈ (0.4, 0.7) where lim m→0 ψ ψ ∞ is plausibly non-zero, implying broken U(2N) symmetry. We therefore deduce the existence of U(2N) symmetry breaking for N = 0 for sufficiently strong coupling, though the data is not of sufficient quality to determine whether there is a critical g 2 c such that symmetry is restored for g 2 < g 2 c , such as occurs in quenched QED 4 with staggered fermions [27]. Nonetheless, this exercise suggests that N c > 0 for the Thirring model.
To emphasise the importance of first taking the L s → ∞ limit, in Fig. 7 we plot the g −2 = 0.4 data vs. m at fixed L s . While the increasing curvature of the data with L s is suggestive, there is no compelling evidence to support a non-zero intercept on the vertical axis as m → 0. The m → 0 and L s → ∞ limits do not commute. For completeness, in Fig. 8 we plot bilinear condensate data for the quenched surface model. Data have been taken at couplings ranging from the relatively weak g −2 = 0.6 to g −2 = 0.2 on 16 3 × L s , with L s = 8 at all couplings, increasing up to 24 (g −2 = 0.3) and 40 (g −2 = 0.2), using at least 12500 HMC trajectories in all cases. The abcissae of some datapoints at this strongest coupling have been slightly displaced for clarity. The emerging picture is qualitatively different from the bulk case; here there is no evidence for any systematic change in the signal as L s is increased. For g −2 ≥ 0.3 the data show a fairly weak g 2 -dependence with lim m→0 ψ ψ = 0 consistent with the absence of symmetry breaking. By contrast data at g −2 = 0.2 admit a plausible extrapolation to a non-zero intercept at m = 0, implying condensation at this strongest coupling, and suggesting both that N c > 0, and that there exists a critical g −2 c > 0 at which the symmetry is restored. Disentangling these effects from the strong-coupling artifacts discussed below Fig. 1 would require an extensive programme of further simulations. To summarize, the quenched exercise shows the importance, at least for the bulk formulation, of taking data at varying L s and performing a plausible extrapolation to the U(2N)-symmetric limit L s → ∞ using the exponential Ansatz (21). While the quenched limit does not correspond to a unitary field theory (the m-independence of the curves in Fig. 6 may be a symptom of this), these results are encouraging because they illustrate at least the possibility of finding symmetry breaking in the Thirring model with DWF. We will apply the same procedure to N = 1 in the next Section.

Results for N = 1
Next we present results from simulations of the full field theory with N = 1, with emphasis on the L s → ∞ limit. As outlined in Sec. 2, the required RHMC simulations are numerically demanding, so the study has been limited to four couplings g −2 = 0.3, 0.4, 0.5, 0.6 chosen to span the region of greatest variation in Fig. 1 while remaining on the weak-coupling side of the maximum. Data were taken at each of 5 masses m = 0.01, . . . , 0.05, and unless stated on a 12 3 spacetime lattice. To probe the large-L s limit we examined L s = 8, 16, 24, 32 and 40, except at the strongest coupling g −2 = 0.3, where runs with L s = 48 were also performed for the three lightest masses. To explore potential volume effects we also simulated a 16 3 lattice at g −2 = 0.3, 0.6 with m = 0.01. For each parameter set a minimum of 600 RHMC trajectories of mean length 1.0 were generated, with observables calculated every five trajectories.

Bilinear Condensate ψ ψ
Just as in the quenched case, the importance of finite-L s corrections increases markedly as the coupling gets stronger. Fig. 9 compares data taken at the strongest and weakest couplings explored for all five mass values, and bears a striking resemblance to Fig. 5. Indeed, the Ansatz (21) again gives a very good description of the condensate data, with χ 2 per degree of freedom for each fit usually < ∼ 2 across the entire dataset. In what follows the condensate values extrapolated to L s → ∞ are based on all available data, with no points excluded from the fit. Fig. 9 also includes results taken on 16 3 denoted by open symbols. On the scale of the plot, volume effects are only discernible at strong couplings and the largest available L s ; fits from both 12 3 and 16 3 will be presented below.  16 3 , and dashed lines fits to (21). The inset plots (ψψ) 2 − ψ ψ 2 vs. L s for various g −2 ; the colour code is that of Fig. 10.
The inset of Fig. 9 shows the variance of ψ ψ , or in physical terms the disconnected contribution to the longitudinal susceptibility, as a function of L s . This demonstrates that the L s → ∞ limit is also key to characterising the fluctuations of a would-be order parameter; indeed it is clear that still larger L s will be needed before this observable converges, particularly at stronger couplings. Note that data from 16 3 are compatible with 12 3 , demonstrating that the observed growth is a finite-L s artifact and is not associated with critical fluctuations. Fig. 10 shows the bilinear condensate ψ ψ following the L s → ∞ limit obtained using (21). Different colours correspond to different couplings -note that uncertainties in the L s → ∞ extrapolation occasionally result in very large errorbars at the weakest coupling g −2 = 0.6. For g −2 ≥ 0.4, the data are consistent with the behaviour ψ ψ(m) ∝ m, implying no symmetry breaking in the limit m → 0. The 16 3 g −2 = 0.6 point suggests that volume effects are small in this regime. This is similar to the findings of simulations with N = 2 using the HMC algorithm [17]; however in that study there was no attempt to take the L s → ∞ limit. Here we rectify that omission by plotting extrapolated data from HMC simulations with L s = 8, . . . , 40 at the strongest available coupling g −2 = 0.3. Fortunately the conclusions of [17] remain unchanged; there is no evidence for spontaneous symmetry breaking, implying N c < 2. The contrast with the results of Fig. 6 is particularly striking; whilst the condensate in the quenched model shows no significant m-dependence, here the linear behaviour is precisely that expected of a unitary field theory in its symmetric phase. For N = 1 at the strongest coupling examined g −2 = 0.3, ψ ψ(m) is a factor of two or greater than data from the next strongest coupling, and a linear extrapolation lim m→0 ψ ψ(m) = Σ ≈ O(0.1) = 0 looks reasonable, particularly if the 16 3 point is used at m = 0.01. This would be consistent with the spontaneous breakdown of U(2) symmetry due to bilinear condensation at this coupling, although non-linear extrapolations to a symmetric limit ψ ψ = 0 cannot at this stage be excluded. If symmetry is indeed broken, on general grounds significant finite volume corrections are expected in the mesoscopic regime mΣV < ∼ 1, and the data support this; note that the dimensionless combination mΣV ≈ 1.5 for the 12 3 , m = 0.01 point.
In summary, Fig. 10 presents strong evidence for the Thirring model with g −2 = 0.3 to exhibit qualitatively very different behaviour from that observed at weaker couplings, due to a significant enhancement of fermion -antifermion pairing. Finite-L s corrections are also much more important in this regime, as illustrated in Fig. 9, and an L s → ∞ extrapolation proves key to interpreting the data. The simplest explanation is that U(2) symmetry is spontaneously broken at the strongest coupling examined, implying N c > 1.

The Approach to L s → ∞
It is interesting to compare the m-dependence of the decay constant ∆, implicitly defined in (21), between different couplings. Of course, for a fixed window in L s , ∆ is easier to pin down for data with large curvature, corresponding to strong couplings and larger masses. For this reason the large uncertainties on ∆ from the weaker couplings g −2 = 0.5, 0.6 don't yield much of use; however results from the stronger couplings g −2 = 0.3, 0.4 plotted in Fig. 11 show a marked contrast. Within sizeable uncertainties ∆(g −2 = 0.4) ≈ 0.06-0.07 is approximately m-independent, whereas ∆(g −2 = 0.3) ∝ m, the linearity becoming more convincing still if the 16 3 value is taken at m = 0.01. The straight line fit shown yields a slope 1.33 (15), with intercept consistent with passing through the origin. This is another hint of a qualitative difference in the behaviour of the model at these two couplings.
Another measure for the approach to the L s → ∞ limit is the residual δ h defined in (14). As shown in [15], it quantifies the difference between the U(2)-equivalent condensates ψ ψ and the measured i ψ γ 3 ψ , and should therefore vanish in a simulation respecting U(2) symmetry. Results for δ h (L s ) for various couplings are shown on a log scale in Fig. 12. Just as in quenched QED 3 (see Fig. 2 of [15]), δ h is strongly couplingdependent. In all cases the data is consistent with an asymptotic behaviour δ h ∝ e −cLs implying U(2) restoration in the large-L s limit; however the restoration becomes slower as coupling increases. There is a marked difference between g −2 = 0.6, where δ h is roughly m-independent, and g −2 = 0.3 where data from all 5 masses are plotted, and c found apparently to decrease systematically with m. At this strong coupling for m = 0.01 δ h is of the same order of magnitude as the signal i ψ γ 3 ψ even for L s = 48. For the larger 16 3 lattice, c is smaller still; a similar trend was observed in [15].
The findings of both Figs. 11,12 are consistent with the extrapolation L s → ∞ used to obtain Fig. 10, and moreover both display qualitative differences between strong and weak coupling, thus supporting the argument that g −2 = 0.3 and g −2 = 0.6 lie in different phases. However the approach to the large-L s limit becomes very slow in the symmetry broken phase in the limit m → 0, which will almost certainly present practical difficulties in future more refined simulations, and may also raise more conceptual problems related to the existence of a U(2)-symmetric limit at strong coupling.

Meson Correlators
Finally we consider correlators of states formed from a fermion and an antifermion, which by analogy with QCD will be referred to as mesons. We will focus on the sector with angular momentum J = 0, in which case four interpolating operators can be written. With a choice of mass term S 3 , they split into two scalars (ψψ,ψγ 3 ψ) and two pseudoscalars (ψγ 5 ψ,ψγ 3 γ 5 ψ). In the event that a symmetry-breaking condensate i ψ γ 3 ψ = 0 forms, then a Goldstone boson of either parity, interpolated byψγ 5 ψ (0 − ) andψψ (0 + ), is expected. The other two states remain massive. The DWF formulation was set out in [15] in terms of "primitive" propagators where the 2+1+1d fermion propagator S(m a ; x, s; y, s ′ ) = Ψ(x, s)Ψ(y, s ′ ) ma . By construction C ±− are real and positive. It can be shown that the meson interpolated bȳ ψγ 5 ψ has a propagator C +− + C −− , while that interpolated byψγ 3 γ 5 ψ has propagator C +− − C −− . The other two mesons in principle require additional fermion propagator calculations with the flip m 3 → −m 3 ; however in the context of quenched QED 3 it was shown in [15] that the propagator interpolated byψψ becomes approximately equal tō ψγ 5 ψ, and that ofψγ 3 ψ equal toψγ 3 γ 5 ψ, in the limit L s → ∞. Degeneracy of these opposite parity mesons is necessary for U(2) symmetry restoration; we will not pursue this issue further, but rather confine our attention to the pseudoscalar channels interpolated byψγ 5 ψ andψγ 3 γ 5 ψ, which will be referred to as Goldstone (G) and non-Goldstone (NG) respectively. Meson timeslice correlators C(τ ) = x C( x, τ ) were calculated on a 12 3 lattice with L s = 40 with m = 0.01 at each coupling already investigated, with a minimum of 500 RHMC trajectories. The primitive propagators were calculated every 5 trajectories by averaging over 5 point sources located at random spacetime points, and the reconstructed G and NG correlators are plotted in Figs. 13 and 14 respectively. Additional calculations with m = 0.05 were performed at g −2 = 0.3, 0.6.
A 12 3 lattice is too far from both thermodynamic and zero-temperature limits for any statements about the model's spectrum to be reliable; there is no sign of pure exponential decay corresponding to a simple propagator pole, and as we shall see below it is also not safe at this stage to infer anything regarding the residue. Accordingly we restrict ourselves to qualitative comments. In the G channel (Fig. 13) there is a huge variation in signal size as the coupling increases, with in particular a factor of 18 increase at the midpoint τ = 6 between g −2 = 0.5 and 0.3, and one of 5 between g −2 = 0.4 and 0.3. Moreover, the impact of changing m from 0.01 to 0.05 is far more pronounced at g −2 = 0.3, where C G (6) decreases by almost a half, and g −2 = 0.6, where the decrease is less than 20%. In the NG channel (Fig. 14) the variation at the midpoint between g −2 = 0.3, 0.6 is an order of magnitude smaller.  Quantitively, the ratio C G (6) : C N G (6) increases from ∼ 1.5 at g −2 = 0.6 to ∼ 17 at the strongest coupling. In terms of primitive correlators, this implies that C −− ≪ C +− at weak coupling, so that G and NG channels are approximately degenerate, but that C −− < ∼ C +− by g −2 = 0.3. This development is also reflected by the relatively large errorbars in C N G at this coupling.
In view of the results for the would-be order parameter ψ ψ of Sec. 4.1, a natural interpretation of these results is that U(2) symmetry spontaneously breaks somewhere in the range g −2 ∈ (0.5, 0.3) and that both the increased magnitude of C G and its enhanced sensitivity to a change in m is due to its developing into a true Goldstone boson.
It is interesting to compare these results with those presented in Fig. 15 of Ref. [17] for the would-be Goldstone meson on 12 2 × 24 in the surface formulation of the Thirring model with N = 2. In that case C G (τ ) manifests a τ -independent plateau for 5 < ∼ τ < ∼ 20 over a range of couplings, interpreted in [17] as being due to fermion propagators reconnecting only after one of them loops around the timelike extent of the system. In other words, the surface formulation does not appear to support mesonic bound states. With the caveats already discussed, the meson correlators of Figs. 13,14 do appear to resemble those of conventional mesons. This is the first hint that the bulk formulation is the preferred approach to the Thirring model with DWF. Finally, we consider more a formulational issue by examining the axial Ward identity, first considered in this context in [17]. For a system with the U(2) symmetry anticipated in the L s → ∞ limit, the following identity relating the order parameter with the integrated correlator holds 1 : The ratio ψ ψ /mχ π is plotted vs. g −2 in Fig. 15, together with bulk formulation results for N = 2 on 12 3 × 16 [17]. The clear issues are that the ratio is neither constant, nor equal to unity, as required by (24). Ref. [17] suggested this is due to a non-trivial relation between either or both of the fermion mass m and the physical fields ψψ defined in (10), and their continuum counterparts (amusingly, on the vertical scale used in Fig. 15 the data provoking this speculation now looks rather constant as a function of g −2 ). This would mean that bilinear operators and/or the fermion mass m would need to be renormalised for the Ward identity to apply. Fig. 15 suggests these considerations become still more important at strong coupling for N = 1; indeed, at g −2 = 0.3 even the effect of changing m results in a marked renormalisation. Recall the ratio was observed to be m-independent for N = 2 [17]. Strong renormalisations depending on both g 2 and m cannot be ruled out for the parameter regime studied in this paper; moreover we draw some encouragement from the hints in Fig. 15 that the effect is smooth as g −2 ranges from 0.6 to 0.3, consistent with UV physics, and in contrast with the sharp changes over the same range reported in the rest of this section, associated with a symmetry breaking phase transition.

Discussion
The main result of this paper is that the spontaneous breakdown of the U(2N) symmetry present for massless reducible fermions in 2+1d can be demonstrated in simulations of an interacting field theory using domain wall fermions. The proof of concept was given in the quenched limit in Sec. 3, where the importance of taking the L s → ∞ limit before the m → 0 limit was shown. Next, simulations of the unitary N = 1 model with a newly-developed RHMC algorithm, discussed in Sec. 4, yielded results following the same procedure consistent with unbroken U(2) symmetry for g −2 ≥ 0.4, but with enhanced bilinear condensation at the strongest available coupling g −2 = 0.3, consistent with a non-vanishing intercept in the m → 0 limit signalling the breaking of U(2) (see Fig. 10). Meson correlators on admittedly small spacetime volumes were consistent with the Goldstone spectrum expected for the breaking pattern U(2)→U(1)⊗U(1) (see Figs. 13,14). The most natural conclusion is that there is a symmetry-breaking phase transition at some g −2 c ∈ (0.3, 0.4), that the critical flavor number in the 2+1d Thirring model satisfies 1 < N c < 2, and that there is the potential for a QCP in the N = 1 model described by a strongly-interacting local unitary quantum field theory. Final confirmation of this important result must await further simulations permitting enhanced control over both V → ∞ and m → 0 limits; until then strictly the bound we have found is 0 < N c < 2. The large disparity with the staggered Thirring model result N c = 6.6(1) [11] is a dramatic indicator of the importance of the faithful rendition of global symmetries when modelling strongly-interacting systems.
As a bonus, the form of the meson correlators strongly suggests the preferred formulation of the Thirring model with DWF uses the bulk formulation of the vector auxiliary, clearing up an outstanding issue from previous work [17]. However the quenched results of Sec. 3 revealing symmetry breaking in the surface model at very strong coupling mandate further investigation of this formulation. The question of the most natural formulation of this (or any strongly-interacting) model permitting systematic numerical investigation remains open; here it is prudent to recall that simulations of the Thirring model with SLAC fermions find N c < 1 [18].
The N = 1 bulk simulations also raise some concerns. The decay constant ∆ governing the L s → ∞ extrapolation seems to follow ∆ ∝ m for g −2 < g −2 c (see Fig. 11), implying that there may be both practical and even conceptual difficulties reaching the massless limit in the broken phase. The residual δ h parametrising the explicit U(2N) symmetry breaking at finite L s is also rather large and slowly-decaying in this regime (see Fig. 12). This suggests rather careful attention will need to be paid to the question of U(2N) symmetry in future studies of the broken phase. We also remark that a further outstanding issue is the locality of the associated 2+1d overlap operator, which governs the validity of the U(2N) in terms of global symmetry rotations on local fields in 2+1d [16,28]. An understanding of each of these issues is a precondition for a satisfactory operational definition of quantum field theories of strongly-interacting fermions.
In future work we plan to implement simulation code with improved performance to counter the considerable numerical effort required for the inversion of M † M at the strong couplings relevant for symmetry breaking. The need for further improvements in the invertor algorthm, and even in the DWF formulation following the ideas of [29], should not be ruled out. The next step is a more refined scan of the N = 1 theory in the critical region g −2 ∈ (0.3, 0.4) with the goal of first locating and then characterising the critical point at g −2 c . The potential difficulty of correctly capturing critical fluctuations is highlighted in the inset of Fig. 9. Finally, using the control over global symmetries furnished by DWF it will be straightforward to examine the effect of a U(2N) and parityinvariant "Haldane" contact interaction (ψγ 3 γ 5 ψ) 2 , which in Ref. [6] was found to be a component of the interaction at the fixed point. Exploratory results in this direction were reported in [26].
In the Dirac basis γ µ = σ µ+1 ⊗ τ 3 (µ = 0, 1, 2) and γ 3 = 1 1 ⊗ τ 2 , where σ and τ are Pauli matrices, and setting m a = 0, we find where B = ∂ 3 , C = (DŴ −∂ 2 3 )∂ −1 3 , and the last step follows if B is invertible. Then detM would be positive definite, and moreover M could be represented as a positive operator making it possible to simulate using bosonic pseudofermions. Now let's examine the commutator. The contributions [∂ 3 , D µ ] = [∂ 3 ,D 2 + M] = 0, which follows provided the link connections obey U µ,x = U µ,x±3 and U † µ,x = U † µ,x±3 . This is the case both for gauge theories and for the bulk formulation of the Thirring model; in each case the connection is "3-static", ie. ∂ 3 U µ,x = 0. However, the remaining part is non-vanishing: Whilst a simple physical interpretation of this term is obscure, it is clear the obstruction to proving the positivity of detM has its origin in the open boundary conditions imposed at the walls. Now consider a Dirac basis γ µ = σ µ+1 ⊗ τ 2 , γ 3 = 1 1 ⊗ τ 3 so that .
In this case the obstruction to proving positivity turns out to be the non-vanishing commutator which by construction is evenly distributed throughout the bulk. The commutator (32) vanishes for configurations in which both the "plaquette" U µν = 1 and U † µ U µ = 1, which for the Thirring model is expected to be reached only in the limit g 2 → 0.
We conclude that detM is real but not in general positive, motivating the use of the RHMC algorithm to simulate the functional measure det(M † M) 1 2 outlined in Sec. 2. We note that there is no such obstruction for HMC simulations of detM using twistedmass Wilson fermions [30] or overlap fermions [31], both of which have been recently used to study QED 3 .