Towards Critical Physics in 2+1d with U(2N)-Invariant Fermions

Interacting theories of N relativistic fermion flavors in reducible spinor representations in 2+1 spacetime dimensions are formulated on a lattice using domain wall fermions (DWF), for which a U(2N) global symmetry is recovered in the limit that the wall separation $L_s$ is made large. The Gross-Neveu (GN) model is studied in the large-N limit and an exponential acceleration of convergence to the large-$L_s$ limit is demonstrated if the usual parity-invariant mass $m\bar\psi\psi$ is replaced by the U(2N)-equivalent $im_3\bar\psi\gamma_3\psi$. The GN model and two lattice variants of the Thirring model are simulated for N = 2 using a hybrid Monte Carlo algorithm, and studies made of the symmetry-breaking bilinear condensate and its associated susceptibility, the axial Ward identity, and the mass spectrum of both fermion and meson excitations. Comparisons are made with existing results obtained using staggered fermions. For the GN model a symmetry-breaking phase transition is observed, the Ward identity is recovered, and the spectrum found to be consistent with large-N expectations. There appears to be no obstruction to the study of critical UV fixed-point physics using DWF. For the Thirring model the Ward identity is not recovered, the spectroscopy measurements are inconclusive, and no symmetry breaking is observed all the way up to the effective strong coupling limit. This is consistent with a critical Thirring flavor number $N_c<2$, contradicting earlier staggered fermion results.


Introduction
Relativistic field theories of particles moving in the plane have received recent attention, particularly within the condensed-matter community, because of their potential role in describing the spin-liquid phase of quantum antiferromagnets [1], the pseudogap phase of cuprate superconductors [2], and of course low-energy electronic excitations in graphene [3]. However, they are interesting to study in their own right. Self-interacting theories of fermions are thought to exhibit an unexpectedly rich variety of ultra-violet renormalisation group fixed points [4], each yielding a new interacting continuum theory. One manifestation of the fixed points is the existence of phase transitions separat-ing massless fermions from a phase in which a mass gap is dynamically generated; in condensed matter physics this represents a transition from a metallic to an insulating phase. Different fixed points fall in different universality classes, which depend on both the pattern of symmetry breaking and the number of interacting species N.
Because the fixed points occur at strong coupling, they present a calculational challenge, exemplified by the generic power-counting non-renormalisability of the perturbative expansion in powers of g 2 for spacetime dimensionality d > 2, since for a four-fermi contact interaction [g 2 ] = 2 − d. At this stage it helps to be more concrete by discussing specific examples. The Gross-Neveu (GN) model for interacting fermions in 2+1d is defined by the continuum Lagrangian density where ψ is an N-flavor 4 component spinor field. The bare mass and interaction terms each reduce the global symmetry from U(2N) to U(N)⊗U(N) (see discussion below (2)); in addition there is a discrete Z 2 symmetry 1 ψ → γ 5 ψ,ψ → −ψγ 5 which is broken by the mass term but not the interaction. Whilst a weak-coupling expansion makes no sense as stated in the previous paragraph, it is possible to develop an alternative expansion in powers of 1/N, favouring Feynman diagrams containing closed loops, and suggesting a resummation [5]. At strong coupling ag 2 ≥ ag 2 c ∼ O(1) in the limit m → 0, where a is a UV regulator length scale, it is found that Z 2 symmetry is spontaneously broken by a vacuum bilinear condensate ψ ψ = 0. There is apparently no obstruction to taking a continuum limit a → 0 as g 2 → g 2 c from either phase. In the resummed theory the interaction is no longer pointlike but rather mediated by exchange of a scalar degree of freedom propagating as k −1 in the deep Euclidean region k → ∞; this anomalous scaling cannot correspond to a term in a local Lagrangian. This persists at higher order in 1/N. Critical exponents receive O(1/N) corrections, but always consistent with hyperscaling, a consequence of 1/N-renormalisability [5,6]. The picture suggested by the 1/N expansion is confirmed by numerical simulations, which observe the symmetrybreaking transition and extract critical exponents for N as small as 2 [7]- [10], or even 1 if a honeycomb lattice is used [11].
The Thirring model has the same global symmetries as N-flavor QED 3 . The Lagrangian (2) is invariant under a U(2N) generated by matrices rotating the N flavors among themselves tensored with the 4 Dirac matrices {1 1, γ 3 , γ 5 , iγ 3 γ 5 }. The parity-invariant mass term mψψ is not invariant under γ 3 or γ 5 rotations, so there is an explicit breaking U(2N) →U(N)⊗U(N). Goldstone's theorem implies the spontaneous breaking of this symmetry results in 2N 2 massless bosons, whereas there are none for the Z 2 GN model of the previous paragraph. However, like the GN model the Thirring model has a renormalisable 1/N expansion [12,13], this time with a resummed vector mediating interactions between conserved currents with UV behaviour ∝ k −1 . The resummation is not associated with a phase transition and the expansion can be developed for any g 2 , implying the coupling is marginal. As g 2 is raised the mass M v associated with the small-k behaviour of the vector propagator varies from 2m at weak coupling to M 2 v ∼ O(m 4−d /g 2 ) at strong coupling [13]. As g 2 → ∞ this suggests (2) is a theory of conserved currents interacting via massless vector exchange, in many respects similar to QED 3 .
However, this may not be the end of the story. Dynamical mass generation through spontaneous symmetry breaking does not occur to any order in 1/N, but for sufficiently large g 2 and sufficiently small N there are self-consistent solutions of Schwinger-Dyson equations which do have this property [14]- [16]. In the limit g 2 → ∞ there is a critical N c below which symmetry breaking occurs: for integer N < N c we therefore expect to find fixed points for some finite g 2 c (N). The problem has also been studied using the functional renormalisation group [4,17], implying the existence of at least two distinct fixed points in the space of possible fermionic theories; however, this approach suggests the nature of the fixed-point interaction is more general than the simple "GN" or "Thirring" forms (1,2) discussed so far, and that a more faithful description requires extra interaction terms consistent with the global symmetries in play. The Thirring model has also been studied using lattice simulations for 2 ≤ N ≤ 18 [18]- [20] which confirm that a symmetry-broken phase is indeed present, that N c ≈ 7 [20], and that critical exponents extracted from the equation of state close to the fixed point depend sensitively on N, quite distinct from the behaviour of the GN model. Since none of these properties is accessed in a systematic weak-coupling method, the 2+1d Thirring model may well be the simplest fermionic QFT requiring a computational approach.
Almost all lattice feld theory studies of 2+1d fermions to date have employed the staggered fermion formulation, in which fields are described by N s -flavor single spinor component Grassmann fields χ,χ defined at each site, and relativistic covariance in the long-wavelength limit ensured by including a space-dependent ± sign on each link with the defining property that the product of such factors around any elementary plaquette equals -1. Well-known algebraic transformations show that a conventional Dirac action is recovered as a → 0 expressed in terms of reducible (ie. 4-spinor) fields ψ,ψ, with ψ defined not at a site but rather distributed over the 2 3 sites defining an elementary cube. Hence ψ is interpreted as describing 2N s flavors of 4-component spinor [21]. However, for a = 0 the staggered formulation does not respect the expected continuum U(4N s ) symmetry but rather a remnant U(N s )⊗U(N s ). Within the lattice community it is widely believed that the full global symmetry is recovered in the weakly-coupled a → 0 limit expected for QCD: there is no reason to believe this is also the case at a stronglycoupled fixed point. For instance, with N = 2 the above considerations suggest distinct breaking patterns of Z 2 for the GN model and U(4)→U(2)⊗U(2) for Thirring. Recent simulations performed with N s = 1 staggered fermions using an efficient fermion bag algorithm, however, have found compatible critical exponents for both "GN" [9] and "Thirring" [22] lattice models, suggesting that for this minimal flavor number the two models describe the same fixed point; in other words, extra microscopic interactions forced by the lower symmetry of the staggered action [18] may be pushing both models into the same renormalisation group basin of attraction. This seems a surprising result when the models are formulated using bosonic auxiliary fields as in the following section, which is both natural for developing the 1/N expansion and required for a conventional hybrid Monte Carlo (HMC) algorithm; however when written purely in terms of χ,χ fields distributed over the vertices of elementary cubes, the interactions differ by only one, presumably irrelevant, term [9].
The considerations of the previous paragraph suggest staggered fermions are not adequate to capture faithfully the correct physics of a fixed point with g 2 c = 0. An approach in which the fermions have the correct global U(2N) symmetry built in is strongly indicated. This insight has been shared by the Jena group, who have recently applied a non-local SLAC derivative operator to the Thirring model [23]. In this paper we will apply domain wall fermions (DWF), originally devised for the study of light quarks in QCD [24,25], and initially studied for 2+1d systems in [26]. The key idea is that a fictitious dimension x 3 is introduced along which fermion propagation is governed by the operatorψ∂ 3 γ 3 ψ. The third dimension has a finite extent L s , with open boundaries called domain walls at each end labelled ±. The only terms coupling the walls are either proportional to the current mass m or are interactions of the GN form (1). Under generic conditions there are exponentially-localised zero-mode solutions of the 2+1+1d Dirac equation at each wall which are eigenstates of P ± = 1 2 (1 ± γ 3 ). It is thus plausible that 2+1d operators and Green functions constructed from 2+1+1d fields living on the walls retain the properties of a theory which is invariant under rotations of the form e iαγ 3 . On a lattice, if the kinetic operator is chosen with a Wilson mass M of opposite sign to m, then the doubler modes generically plaguing lattice fermion formulations do not couple to normalisable modes and are hence irrelevant [24]. Moreover, in the large-L s limit it has been shown both numerically [26] and analytically [27] that full U(2N) symmetry is recovered, ie. e iβγ 5 and e −δγ 3 γ 5 rotations also become invariances. An unanticipated bonus is that the approach to the large-L s U(2N)-invariant limit is accelerated if instead of the hermitian mψψ the physically equivalent, but antihermitian, mass term im 3ψ γ 3 ψ is used. This playoff between the different forms of parity-invariant mass term available for reducible spinor representations in 2+1d has also recently been exploited in a lattice study of non-compact QED 3 [28].
The remainder of the paper is organised as follows. In Sec. 2 we review the DWF formulation for 2+1d reducible fermions, setting out the different possible parity-invariant mass terms and the approach of the corresponding bilinear expectation values to the large-L s limit first studied in the context of quenched QED 3 in [26]. Lattice versions of the GN and Thirring models using DWF are then proposed, and their simulation using an HMC algorithm outlined. While the lattice transcription of the GN model is fairly straightforward, based on a bosonic scalar auxiliary field confined to the walls [29], there are (at least) two possible ways to treat the Thirring model, one in which a vector auxiliary field A µ is confined to the walls, and one in which A µ is defined uniformly throughout the bulk 0 ≤ x 3 ≤ L s in analogy with the treatment of gluon degrees of freedom in QCD with DWF.
Next, Sec. 3 examines the GN gap equation in the large-N limit, which predicts a fixed point and spontaneous dynamical mass generation for g 2 > g 2 c . We build on the pioneering work of Ref. [29] by generalising their solution to the case of the antihermitian mass term im 3ψ γ 3 ψ and demonstrating an exponential improvement in convergence to the large-L s limit as a result. The gap equation also serves as a check to simulations of the GN model with N = 2 presented in Sec. 4, where results for the dynamically-generated gap Σ(g 2 ) are used to monitor the approach to the large-L s limit, and equivalence of the hermitian and antihermitian mass terms is demonstrated. The essential question of whether symmetry breaking and critical behaviour can be probed using DWF is also addressed through studies of the scalar susceptibility peaking in the vicnity of the critical point, recovery of the axial Ward identity (using a variant of the GN model with a U(1) axial symmetry), and finally for the first time the fermion propagator obtained with DWF is used in an exploratory study of both fermion and meson masses. While no attempt is made at a full-blown characterisation of the nature of criticality, the results of this section support the physical scenario outlined above and captured in the large-N expansion [6].
In Sec. 5 we turn attention to the Thirring model with N = 2, presenting results of HMC simulations of both surface and bulk models. The symmetry-breaking bilinear condensate is calculated as a function of g 2 and evidence presented that physical couplings all the way up to the strong-coupling limit are probed. Each model shows evidence for the influence of interactions, as does the auxiliary boson action, but the results differ in qualitatively important ways. Significantly, the axial Ward identity is not respected, signalling that the relation of lattice fields and parameters to the putative continuum theory is not yet under control, and that at this stage it is not yet possible to state whether surface or bulk approaches is optimal. Fermion spectroscopy in this case is hindered by large phase fluctuations, whereas meson spectroscopy requires a much larger temporal extent than the L t = 24 studied here. A robust finding, however, is that there is no evidence for spontaneous symmetry breaking, ie. our results support lim m→0 ψ ψ = 0, on volumes and using lattice parameters where symmetry breaking is clearly observed using staggered fermions [18]. This strongly suggests that for a theory of the form (2), N c < 2. The results are summarised and discussed in Sec. 6, and an Appendix contains technical details of the free fermion propagator using DWF, needed for the large-N calculation of Sec. 3.

Formulation and Simulation
First, let's define the lattice action to be studied. The fermion kinetic term uses the 2 + 1d domain wall operator defined in [26,27]: where the fields Ψ,Ψ are four-component spinors defined in 2+1+1 dimensions. The first term D W is the 2 + 1d Wilson operator defined on spacetime volume V with M the domain wall height parameter, and D 3 controls hopping along the third dimension separating the domain walls at s = 1 and s = L s : Here the projectors P ± ≡ 1 2 (1 ± γ 3 ) and the connection link fields U µ will be specified more fully later, with U µ ≡ 1 for free fields. The mass term in (3) only involves fields on the domain walls themselves, and can be chosen as a linear combination of terms which are either hermitian: or antihermitian: or m 5 S 5 = im 5 xΨ (x, 1)γ 5 P − Ψ(x, 1) +Ψ(x, L s )γ 5 P + Ψ(x, L s ).
Next consider 2+1d fields ψ,ψ defined on the walls as follows: in terms of which the three mass terms are written m hψ ψ, im 3ψ γ 3 ψ, and im 5ψ γ 5 ψ respectively. The three terms are all parity-invariant and define physically equivalent ways of breaking U(2N) →U(N)⊗U(N). In Ref. [26] it was demonstrated, in the context of quenched QED 3 , that for sufficiently large L s the three mass terms yield compatible results for the corresponding bilinear condensates ψ Γ i ψ , with Γ i ∈ {1 1, iγ 3 , iγ 5 }, consistent with the recovery of U(2N)-invariance in this limit. Moreover, while i ψ γ 3 ψ and i ψ γ 5 ψ are numerically indistinguishable, the finite-L s errors show a distinct hierarchy, ie. with then ∆ h ≫ ǫ h ≫ ǫ 3 ≡ ǫ 5 . The error ∆ h is defined by the imaginary part of i ψ γ 3 ψ obtained using justΨ(L s ) and Ψ(1): the fields on the opposite walls yield the conjugate. Therefore ∆ h can be estimated using measurements made with mass term m 3 S 3 . In Ref. [27] it was shown for a gauge theory that an expansion of the bilinear condensate ψ ψ Ls in powers of m h /D Ls , where D Ls is a 2+1d truncated overlap operator proportional to the continuum Dirac operator in the large-L s long-wavelength limit, for finite L s generically contains even powers of m h , whereas the corresponding expansion of i ψ γ 3 ψ only contains odd powers of m 3 , a property shared with the continuum theory. Hence a residual ∆ h (m h , Ls) with only weak dependence on m h as m h → 0 cannot be excluded, consistent with the hierarchy reported below (10).
The Gross-Neveu (GN) model for interacting fermions in 2+1d, defined by the continuum Lagrangian density (1) with ψ an N-flavor 4 component spinor field, exhibits spontaneous breaking of Z 2 . The model is readily generalised to exhibit spontaneous breaking of a continuous symmetry, eg. U(1), by modifying the contact interaction to [(ψψ) 2 − (ψγ 5 ψ) 2 ] -see Sec. 4.2 below. It is convenient to reformulate (1) in terms of a real scalar auxiliary boson field σ: in which case symmetry breaking is signalled by Σ ≡ σ = 0. Note that physically equivalent models are obtained by replacing the contact interaction of (1) by the U(2N)equivalent forms −(ψγ 3 ψ) 2 , −(ψγ 5 ψ) 2 , along with masses m 3 , m 5 multiplying the corresponding bilinears.
Formulation of the GN model on a lattice with DWF proceeds from the observation that the interaction with the auxiliary in (11) formally resembles a mass term [29]. The DWF formulation follows from (4) with U µ ≡ 1, (5),(6) with the interaction term defined solely in terms of fields on the domain walls, with obvious generalisations based on (7,8). An interesting distinction with previous work is that here the auxiliary field variables are simply defined on the lattice sites; in the conventional formulation using staggered fermions σ is defined on the sites of the dual lattice [6]. The other interacting theory considered in this paper is the Thirring model (2). Again, it is convenient to recast the model using a real vector auxiliary field A µ : In this latter form the formal resemblance of A µ to an abelian gauge field is manifest, although the last term in (13) spoils gauge invariance. In the 1/N expansion A µ interpolates a massive vector boson of mass M v ; the ratio M v /m is governed by the coupling strength g 2 [13]. However, symmetry breaking U(2N) →U(N)⊗U(N) via generation of a bilinear condensate does not occur at any order in 1/N. There are several variants of lattice formulation of the Thirring model, even when using staggered fermions [18]. In the so-called non-compact approach the interaction between fermion bilinears and the vector auxiliary defined on the lattice links is linear; this has the virtue that only four-fermion terms are generated on integration over A µ , making the connection with the continuum form (2) as transparent as possible. However, as shown in [18], this regularisation fails to preserve transversity of the vacuum polarisation term contributing to the A µ -propagator (ie. ∂ µ Π µν = O(a −1 )), leading to an additive renormalisation of g 2 and consequent uncertainty in identifying the strongcoupling limit [20]. In this paper two non-compact DWF formulations are investigated. First, by analogy with (12) we study a surface formulation with U µ ≡ 1 in (4) and link fields A µ (x) defined solely on the walls interacting with point-split bilinears: Notice in this case the interaction couples fermion fields on the same wall. Second, we push the analogy between the vector auxiliary and an abelian gauge field by defining a bulk interaction between an s-independent A µ (x) and the vector bilinear current defined for all s: (15) This differs from (14) on the walls by the presence of a formally-irrelevant remnant of the Wilson term (corresponding to the ±1s in (15)). The relation with the gauge-invariant kinetic term (4) with U µ = (1 + A µ ) is clear. If we regard A µ as a gauge field, then the distinction between (14) and (15) is that in the former case the 2+1+1d fields are exposed to s-like plaquettes carrying non-zero flux at both s = 1 and s = L s , whereas in the latter case such plaquettes carry zero flux by construction. At strong coupling the analogy may be crude; since the effective connection is (1 + A µ ) rather than e iAµ the s-plaquettes are not constrained by unitarity, and may still fluctuate in magnitude if not in phase.
After introduction of auxiliary bosons both GN and Thirring models with DWF can be written in the form ie. bilinear in the 2+1+1d fields Ψ i ,Ψ i , where explicit flavor indices are shown. For the bulk Thirring model interactions are encoded within S kin and there is no separate S int . The interaction S int is one of (12,14,15) and the bosonic action S bos is an obvious lattice generalisation of the quadratic terms in (11,13). On the assumption that M, M † describe similar physics, then the HMC algorithm may be used to simulate both models starting from the equivalent pseudofermion action The requirement to have a positive definite kernel means that N must be chosen even, and hence the minimal number of flavors simulable with the HMC algorithm is N = 2. However, in order to obtain the correct functional measure, in general it is necessary to correct for the effect of unphysical bulk fermion modes [25], so that the fermion operator coincides with a 2+1d overlap operator in the large-L s limit. For U(2N)invariant fermions this is done by including for each flavor a term det(D −1 (m h a = 1)) in the functional measure [27], which may be thought of as arising from integration over bosonic Pauli-Villars fields with action ζ † D(1)ζ. It is computationally efficient to use the same pseudofermion fields Φ, Φ † for both fermions and Pauli-Villars fields, and the following action results: Since the GN and surface Thirring models are formulated with U µ = 1 in (4), the Pauli-Villars kernel D(1) has no dependence on the bosonic variables, and so can be dropped with no dynamical impact. Hence in these cases the simpler form (17) may be used, as pointed out in [29]. For the bulk Thirring model, the action (18) is simulated. For all the results presented in this paper the domain wall height is chosen to be Ma = 1.

Insights from Large N
One of the main questions to be addressed in this paper is how critical physics appears when DWF are used, and what is the resulting dependence on additional non-physical parameters introduced in the formulation such as L s . The GN model provides a good starting point because critical behaviour is already manifest in the large-N approximation, and can be accessed analytically. This approach was first applied using DWF in [29]; the corresponding study for staggered fermions was performed in [6]. We start from the continuum GN model defined in (11). Spontaneous breaking of a Z 2 global symmetry is signalled by the development of a vacuum expectation Σ = σ , which in the large-N limit is given self-consistently by the gap equation For DWF σ is localised on the walls according to (12), and the gap equation becomes the subscript h denotes that we initially focus on a hermitian interaction term, with the fermion propagator given by D † G, and the free fermion Green function G(p; s, s ′ ) where p is a 2+1d momentum is derived in Appendix A. Throughout this section units are defined such that a = 1. The first term in square brackets contains the product D † (1, s)G(s, L s ). Using (51,52) we find The coefficients B, A ± , A m are given in (55) In calculating P − (D † G)(1, L s ) and P + (D † G)(L s , 1), terms proportional top / can be dropped since they vanish on tracing. The resulting gap equation is where the mode sum on an L x L y L t lattice runs over p x,y = 2πn x,y /L x,y , p 0 = 2π(n 0 + 1 2 )/L t . Note that finite-L s corrections appear at both O(e −α(Ls−1) ) and O(e −2α(Ls−1) ).
In the limit L s → ∞, we take first the massless limit m h → 0 and then the limit Σ → 0 to find the critical coupling, using (57): The summand is given, using (22), by where the factor z(p) = 1 − be −α was introduced in [29].   = 0.408523 from a massless phase to a phase with spontaneous mass generation given by the solution of (25). It can be seen that the finite-L s corrections are significant for L s < ∼ 6, and still discernable even for L s = 12; note that Σ h (g 2 ) approaches the large-L s limit from above.
With mass term proportional to m 3 the gap equation becomes NΣ 3 /2g 2 −i ψ γ 3 ψ = 0 and (21,23) are replaced by and The result is with with A m 3 given by (66); the full gap equation now reads In the large-L s limit (32) coincides with (24,25) as it must; remarkably, however, since B is free of finite-L s corrections, and from (67)  Solutions of (32) are also plotted in Fig. 1. As predicted, the finite-L s corrections are much smaller, and essentially under control by L s = 6. Also note Σ 3 (g 2 ) approaches the large-L s limit from below. This corroborates the improved properties of the "twisted mass" formulation with respect to approaching the U(2N)-symmetric limit at large L s , observed empirically in quenched QED 3 in [26], and demonstrated analytically for gauge theories in [27]. Finally, Fig. 2 shows the approach to the large volume limit for fixed L s ; as the volume increases the expected scaling Σ 3 ∝ (g −2 c − g −2 ) is recovered in the symmetry-broken phase, consistent with the large-N critical exponent β = (d − 2) −1 + O(1/N 2 ) [6]. As expected, finite volume effects become significant for Σ 3 < ∼ L −1 x .

Numerical Results for the Gross-Neveu Model
For finite N the results of the previous section are subject to quantum corrections. In principle for the GN model these are calculable via the 1/N expansion, but we will use numerical simulations to address the question of what critical behaviour of DWF fermions looks like under these circumstances. The results presented in this section were obtained using a HMC algorithm based on the action (17)   Initially we focus on the GN model with Z 2 symmetry defined by the continuum action (11). The order parameter Σ = σ is related to the corresponding single-flavor bilinear condensate by an equation of motion Σ h = g 2 ψ ψ . The bilinear condensate may be estimated by stochastic means [26] at a significant numerical overhead. We have checked that the simulation respects the equation of motion, but chose to focus resources on the observable Σ. Fig. 3 shows results for both Σ h and Σ 3 as functions of g −2 on a 12 3 × L s system, with bare mass am = 0.01. For L s = 2 Σ h is in approximate agreement with the large-N result of Section 3, but as L s increases the trend is for Σ(g 2 ) to fall below the large-N prediction as a result of quantum corrections. As before, the finite-L s artifacts for Σ 3 are clearly smaller than those of Σ h , but in this case Σ 3 approaches the large-L s limit from above. There is fairly rapid convergence to the large-L s limit: Σ h (L s = 8) ≈ Σ 3 (L s = 4). This can be checked in closer detail in the right panel. The Σ 3 (L s = 20) data, shown on both panels, were obtained from 20 -30×10 3 HMC trajectories of mean length 1.0 and can be taken to define the effective large-L s limit. The L s = 12 results for both Σ h and Σ 3 are consistent within statistical errors.
The pronounced kink in Σ 3 (L s = 20) seen in the right hand panel of Fig. 3 hints at a critical point at ag −2 c ≈ 0.32 -0.34, well below the large-N value predicted by (25) as expected. Of course, a good estimate of g −2 c requires demonstrable control over finite volume artifacts and the m → 0 extrapolation, but such a simulation campaign is beyond the scope of the current study. We can note, however, that finite-N corrections are large, being O(50%) in the critical region.
The slight increase in the size of the errorbars in the critical region just discernable in Fig. 3 is a signal of critical fluctuations, which are more properly quantified by the susceptibility χ = (σ − σ ) 2 , plotted in Fig. 4. The peak in the critical region signals divergence in the infinite volume massless limit, where we expect χ ∝ |g −2 c − g −2 | −γ with γ = 1 + O(1/N) [6]. The contrast is striking: χ h shows no sign of critical behaviour for L s ≤ 4, and approaches the large-L s limit from below, whereas χ 3 approaches the limit from above. The smallest value where the two are plausibly consistent is L s = 12, but even for L s = 20 there are small differences in the data. Unfortunately, the calculation is hard to control, particularly on the strong-coupling side of the transition, due to occasional brief tunnelling between true and false vacuum states related by Z 2 symmetry, which in this context should be regarded as a finite volume artifact. For this reason the actual peak height in the large-L s limit is hard to estimate from Fig. 4. It appears that near criticality the susceptibility presents a more stringent challenge to reaching the large-L s limit than the order parameter.

Axial Ward Identity
Continuous global symmetries in field theories imply the existence of Ward identities relating Green functions. If we wish to check restoration of a symmetry which is formally broken at the Lagrangian level, it is important to examine the recovery of Ward identities, to check both the symmetry itself and the applicability of field identifications such as (9). To follow this agenda in the GN model it is necessary to enhance the model by changing the broken symmetry from Z 2 to U(1), requiring the introduction of a second bosonic pseudoscalar auxiliary field π. The continuum Lagrangian becomes and the U(1) symmetry We will focus on the axial Ward identity ψ ψ where all Green functions are normalised to just a single fermion flavor. In a symmetry broken phase with lim m→0 ψ ψ = 0, the resulting divergence of the RHS signifies the Goldstone nature of the π field. In the GN model, the Goldstone mode is dominated by disconnected fermion-line diagrams [30]. However, the auxiliary equations of motion may be used to recast the identity as where χ π is the transverse susceptibility (π − π ) 2 , so all required expectation values involve solely bosonic fields. Finally, note that if instead the mass term im 3ψ γ 3 ψ is chosen then the interaction in (33) takes the form iψ(γ 3 σ − γ 5 π)ψ, but the same Ward identity (36) results. Fig. 5 shows data from the N = 2 GN model with U(1) symmetry and mass term m 3 S 3 on a 12 2 × 24 spacetime lattice with L s = 8. The data results from 30000 HMC trajectories of mean length 1.0. Each side of Eqn. (36) is plotted separately; on this scale the errors in Σ are hard to discern whereas the χ π data suffer from large fluctuations due to the Goldstone nature of π, similar to the staggered fermion observations of [30]. Nonetheless, within the admittedly large errors the data are consistent with the Ward identity (36).

Spectroscopy
Finally we present some exploratory spectroscopy. Since the GN model is not constrained by Elitzur's theorem, it is possible to study the propagator of a single fermion. In addition we will examine the simplest meson correlator formed from connected fermion lines, which as shown in [26] interpolates states with J P = 0 ± . This study uses the same ensembles as Sec. 4.2.
First consider the timeslice propagator of a free fermion with mass m f : with P 0± ≡ 1 2 (1 ± γ 0 ) and the sign chosen according to the sign of the temporal displacement x 0 . Using the identification (9) and the identities P ± P 0+ P ± = 1 2 P ± , P ∓ P 0+ P ± = 1 2 γ 0 P ± , the corresponding 2+1+1d correlator with mass m h and x 0 > 0 is trP 0+ (P − Ψ(0, 1) + P + Ψ(0, L s ))(Ψ(x, L s )P − +Ψ(x, 1)P + ) = (38) The generalisation to mass m 3 is straightforward. We have measured timeslice correlators using the first and third terms of the RHS of (38) using 5 randomly-located sources on configurations separated by 5 HMC trajectories (taking care to correct for anti-periodic temporal boundary conditions when x 0sink < x 0source ). They yield two distinct estimates of the fermion correlator labelled 1 1 (formed from 2+1+1d propagators linking the two domain walls) and Γ 0 (formed from propagators starting and ending on the same wall) in the following. The two should coincide in magnitude if the correlator is dominated by a simple pole of the continuum form (37).   6 shows the raw correlators C Γ 0 (x 0 ) and C 1 1 (x 0 ) obtained at coupling ag −2 = 0.24 using both a point source and a Gaussian smeared source with D ⊥ the spatial part of (4). We chose c = 0.25 and N smear = 10. The essential feature is that C Γ 0 is even about x 0 = L t /2, whereas C 1 1 is odd, so that their linear combination is not symmetric about the centre of the lattice. Lines are drawn though the point source data to emphasise this, which also suggest the two channels don't yield signals of equivalent magnitude. The data obtained with smeared sources has the same symmetry, but with much larger errors. These originate in the large phase fluctuations of the background auxiliary Φ field, and render the raw correlators useless for the precision fitting required by spectroscopy with the available statistics. The fluctuations also afflict the point-source data; though invisible on the scale of Fig. 6, the C 1 1 datapoints actually have the wrong sign near the centre of the lattice.
The pragmatic solution adopted here is instead to study the functionsC i = C * i C i , effectively ignoring the phase fluctuations. It should be borne in mind thatC thus defined is not a Green function, and that any resulting particle mass estimate must strictly be a lower bound. It is also worth noting that mass fits to fermion correlators .36 One s Thirring g -2 =0.8 in the U(1) GN model with staggered fermions were obtained without the need for this step [31]. Results forC Γ 0 andC 1 1 for two representative couplings are plotted on a logarithmic scale in Fig. 7, and clearly suffer far less from fluctuations. By constructionC is symmetric about the lattice midpoint. The correlators evaluated with point sources show no coupling dependence for |x 0 | < ∼ 5; we ascribe this to the influence at short temporal separations of excited states which are probably lattice artifacts. This is in notable contrast to using staggered fermions, where excited states are absent permitting fitting over almost the entire temporal extent, eg. [18]. For this reason fits been made using smeared sources, which yield correlators with a better projection onto the ground state and showing a much cleaner g −2 -dependence; reasonable fits to a simple pole were found for x 0 ∈ [6, 18] (Γ 0 ) and x 0 ∈ [5,19] (1 1), and two such fits toC Γ 0 are shown. Most of the data of Fig. 7 were obtained with L s = 8 and mass term m h S h ; we also show one correlator using L s = 16 and m 3 S 3 , and it is clear that at this level of accuracy the large-L s limit is secure. The smeared source data in Fig. 7 also show thatC 1 1 ≈C Γ 0 over the whole x 0 range, but that there are systematic differences near the lattice midpoint, which may be aC-artifact; C 1 1 should vanish at the lattice midpoint, and is ideally fitted with an odd function. Fig. 8 shows results for the fermion mass m f for ag −2 ∈ [0.24, 0.36] and am = 0.01, .02, .03, together with the Σ data from Fig. 5. In the large-N limit m f = Σ + m [6]. The plot shows that this relation is approximately observed, but the measured m f falls systematically below Σ at strong coupling and above Σ at weaker coupling. The fits also yield m f Γ 0 < m f 1 1 , with the trend becoming more marked at weak coupling, as might be anticipated from Fig. 7. Both m f and Σ show similar variation with m over the whole range studied. Determining whether the origin of the mismatch is due to finite spatial volume, the fitted x 0 -range, O(1/N) corrections, or an artifact of fittingC rather than the Green function C is beyond the scope of this exploratory study. Fig. 8 also shows mass fits to meson correlation functions, defined by the combination C +− + C −− using the notation of Sec. 5.2 of Ref. [26] (see (42) below), corresponding to states interpolated by the bilinearsψγ 5 ψ (J P = 0 − ), andψγ 3 ψ (0 + ). Again, there is evidence of significant excited state contamination (see Fig. 15 below), and the fits shown here were obtained from x 0 ∈ [6,18]. There was some difficulty in obtaining stable fits at strong coupling using smeared sources, but by and large point and smeared sources yield compatible results. All the previous remarks about systematic effects apply here; the main feature revealed in Fig. 8 is M 0 ∓ ≈ 2m f . Although the Goldstone mode has quantum numbers 0 − , it is only accessed via disconnected fermion line diagrams [30], or perhaps more effectively via the auxiliary π field as in Sec. 4.2. Mesons formed from connected lines are only weakly bound by O(1/N) effects, hence the spectrum revealed in Fig. 8 is physically plausible.

Numerical Results for the Thirring Model
Next we turn to the Thirring model, for which the bosonic auxiliary field A µ is not simply related to a bilinear condensate order parameter, and where there is no straightforward analytic approach to compare with numerical results. Moreover as discussed in Sec. 2 the lattice prescription is not unique. Accordingly we will explore the models defined by both surface (14) and bulk (15) interaction terms. Unless otherwise stated, the results of this section were obtained with 5000 HMC trajectories over a range of couplings g −2 on a 12 3 system with L s = 16 and am 3 = 0.01, with N = 2 fermion flavors. The residual ∆ h (L s = 16) defined in (10) ranges between 0.6 -1.0×10 −6 so the results are safely in the large-L s limit implying ψ ψ = i ψ γ 3 ψ . For reference, Fig. 9 plots the mean number of congugate gradient iterations to achieve a residual norm of 10 −9 per vector component, needed in the HMC acceptance step, for each model as a function of g −2 .
The relative cost of the bulk model rises steeply as the coupling becomes strong. For . The expectation value is measured using 10 stochastic estimators every 5 HMC trajectories, as described in Sec. 5.1 of [26]. Over the couplings explored the condensate varies by about 20%, and shows marked non-monotonic behaviour, peaking at ag −2 ≈ 0.2. The plot also includes data taken with L s = 20, showing that the large-L s limit is effectively reached, and data taken on 16 3 and 20 3 systems (the latter with 2000 HMC trajectories) showing significant volume effects, though smaller than those shown in Fig. 2 for the large-N GN model in the critical region.
The peak in the order parameter at strong couplings has also been observed in simulations of the Thirring model with staggered fermions [18,20]. In [18] it was observed that the fermion-auxiliary interaction fails to preserve transversity of the vector current correlator, ie. where Π µν is the vacuum polarisation tensor. Transversity originates in Ward-Takahashi (WT) identities arising from an underlying gauge symmetry, which on a lattice implies the link field is represented by e iAµ rather than simply A µ . In QED 3 the WT identity follows from a cancellation of an O(a −1 ) divergence between the two diagrams shown in Fig. 11. With a linearised interaction of the form (14) the right hand diagram is absent because there is no 2-fermion 2-boson vertex. The resulting linear divergence is absorbed by an additive renormalisation of the coupling: g −2 R = g −2 − J(m, N)a −1 . The physical strong coupling limit g −2 R → 0 is thus found at non-zero g −2 ; in practice its location must be determined by numerical simulation [20]. For g −2 R < 0 the vector correlator in the 1/N expansion becomes negative, signalling violation of reflection positivity. Now, the WT identity is independent of the details of the lattice fermion regularisation; even without a detailed calculation of the diagrams in Fig. 11 using DWF it is reasonable to apply the same arguments to the current case. Hence we interpret the peak in Fig. 10 as evidence that the effective strong coupling limit lies at ag −2 ≈ 0.2, and that the simulations have thus explored a range of couplings up to this limit. The variation of i ψ γ 3 ψ with g −2 shows clear evidence for interaction effects. We now observe that data taken with am 3 = 0.005, 0.01, 0.02 lie on top of each other, or in other words, there is no evidence to contradict the hypothesis that lim m 3 →0 i ψ γ 3 ψ = 0 for all values of the coupling. This is in marked contrast with results obtained using staggered fermions on the same volume with comparable lattice parameters; compare Fig. 7 of Ref. [18]. We conclude that a spontaneous symmetry breaking U(2N) →U(N)⊗U(N) is absent in the Thirring model defined by (14) for N = 2.  Fig. 12 shows the results of a similar study for the bulk model (15), this time using the pseudofermion action (18) to perform the HMC simulation. The magnitude of i ψ γ 3 ψ /m 3 is considerably larger, reflecting the fact that the two lattice models are different regularisations of a field theory. Again, there is evidence for g −2 -dependence, and a local maximum at ag 2 ≈ 0.2, this time followed by a steep rise at stronger couplings. Data taken at different m 3 lie on top of each other following rescaling, once again consistent with the absence of symmetry breaking. An interesting contrast between the two formulations is highlighted in Fig. 13 plotting the boson action g −2 A 2 µ per lattice site. For non-interacting fields the expected value is 3 2 . In the surface model the action density stored in the auxiliary fields exceeds the free-field value and increases with coupling strength, whereas the bulk model exhibits the opposite trend, starting from the right below the free-field value and decreasing up to the effective strong coupling limit at ag −2 > ∼ 0.2. Large UV artifacts might be expected for the expectation value of a composite operator, and indeed this is the preferred interpretation for what are ostensibly two different regularisations of the same theory. Nonetheless, the contrast between surface and bulk models may prove a useful diagnostic. For comparison the corresponding quantity g −2 σ 2 is plotted for the Z 2 GN model of Sec. 4. Here there is a clear distinction between near free-field behaviour at weak coupling and a sharp upward rise in the symmetry-broken phase, readily understood since σ is also an order parameter field.
Next consider the axial Ward identity as test of the extent to which U(2N) symmetry is restored. The equivalent identity has been found to hold in simulations of the Thirring model with staggered fermions [19]. For a U(2N)-invariant theory such as the Thirring model in the limit m → 0, the axial Ward identity (35) generalises to where no sum is implied by repeated explicit flavor indices. The mesons interpolated bȳ ψγ 3 ψ,ψγ 5 ψ have opposite parities. With mass term m 3 , the equivalent identity hasψγ 5 ψ as the pseudoscalar, and the equivalent correlator contains contributions from 2+1+1d propagators S(m 3 ; 0, s; x, s ′ ) = Ψ(0, s)Ψ(x, s ′ ) both running between the walls, and starting and ending on the same wall [26]: Assuming that only connected fermion line diagrams contribute to the Ward identity, we define the pion susceptibility χ π = x C 3 π (x). Fig. 14 plots the ratio i ψ γ 3 ψ /m 3 χ π as a function of g −2 , which should take the value unity if the axial Ward identity is preserved by the lattice formulation, for both surface and bulk models. There are clear problems both in terms of the ratio's magnitude and also its variation with g −2 , which is smooth but significant in the regime g −2 R > 0. It is interesting that the trend is opposite for surface and bulk models, again suggestive that the g −2 variation is a UV artifact. There is no variation with m 3 .
The Ward identity is not respected because the bare action (3,14,15) is not U(2N)invariant. Possible causes of the breakdown could be that the correct fermion mass in (41) is not simply related to the lattice parameter m 3 , or that the field identification (9) needs modification, resulting in renormalisation of fermion bilinears, once interactions are present. Whilst these are not fatal objections, they do make it clear that care will be needed in applying DWF techniques to this strongly-interacting system. In particular, Fig. 14 provides little guidance as to whether to choose bulk or surface formulations for further study. One possible way forward is instead to regard the Ward identity as a relation between renormalised quantities, so that m in (41) is replaced by m f , which as a spectral quantity is much better-defined. The physical fermion mass was successfully measured in Thirring model simulations using staggered fermions [18]. To this end the fermion propagator on a 12 3 × 24 lattice with L s = 16, am 3 = 0.01 was studied using 45000 HMC trajectories of the surface model, with measurements made every 5 trajectories using 5 randomly chosen sources. The best results were obtained with a smeared source (39) with D ⊥ incorporating a link connection of the form U µ = e iAµ . The resulting C Γ 0 (ag −2 = 0.8), where positive, is plotted in Fig. 7 showing that in contrast to GN a correct treatment of phase fluctuations is essential to capture Thirring dynamics signal, the fluctuations are still too large to permit a credible fit for m f ; this may well reflect the resemblance of the Thirring model to a gauge theory, for which the correlator vanishes unless a gauge-fixing is specified. Finally, Fig. 15 plots the meson correlator C +− + C −− , obtained with point sources, for the surface Thirring model at two representative couplings. Corresponding results for the GN model discussed in Sec. 4 yielding the spectrum plotted in Fig. 8 are shown for comparison. The contast is clear; the GN data permit a simple pole fit describing a massive meson as reported in Fig. 8, while the x 0 -independent plateaux in the Thirring correlators are due to fermion propagators reconnecting after looping around the timelike extent of the system, and are characteristic of non-confining theories containing light fermions. Extracting spectral information in the meson channel would require lattices of much greater temporal extent than those studied here. Fig. 15 provides further indirect evidence that over the coupling ranges explored dynamical fermion mass generation is happening in the GN model, but not in the Thirring model.

Discussion
Let us summarise the main results of the paper. The programme to apply DWF to relativistic fermions in reducible spinor representations begun in [26,27] has been developed to cover non-perturbative simulations of interesting quantum field theories. The main results of the earlier work, namely that U(2N) global symmetries are recovered in the large-L s limit, and that approach to the large-L s limit is accelerated if an antihermitian or "twisted" mass term im 3ψ γ 3 ψ is chosen, have been confirmed. The large-N solution of the GN model presented in Sec. 3 provides a particularly nice illustration, since here the acceleration is actually exponential. Simulations of the GN model set out in Sec. 4 provide qualitative support for the physical picture revealed by the large-N limit, but also enable a quantification of quantum corrections of O(1/N). Crucially, the results for the gap Σ and scalar susceptibility χ demonstrate that critical physics can be observed using DWF, and that finite-L s artifacts can be controlled. Whilst a quantitative understanding of the critical properties and universal features of the fixed point theory would require simulations on a much larger scale, and in particular require much larger spacetime volumes, there is no reason to doubt the feasibility of such a campaign.
From a theoretical perspective, the principal result of the paper is that with N = 2 the physics of the GN and Thirring models is very different, in contradiction to results obtained with staggered fermions [9,22]. The most obvious distinction is that the GN model exhibits a phase transition at strong-coupling to a phase in which a global symmetry (Z 2 in the example studied) is spontaneously broken and a fermion mass dynamically generated; no such transition is observed in the Thirring model despite strong evidence that the physical strong coupling limit is probed. In the GN case, subsidiary measurements of the axial Ward identity and the mass spectrum yielded results consistent with large-N expectations. This success may be due in part to the fields σ and π being related via equations of motion to bilinears of direct interest such as the order parameter field; it is certainly the case that sampling the {π} ensemble is a very effective means of estimating correlators formed from disconnected diagrams. To our knowledge the measurement of the fermion correlator is the first using DWF; compared to what is known using staggered fermions, the contributions of excited state artifacts are surprisingly large, necessitating an approach based on source smearing. This in turn exacerbates the influence of phase fluctuations, so that the analysis needs to be based on the quasi-Green functionsC. It is remarkable that even so the resulting spectrum shown in Fig. 8 matches large-N expectations as well as it does.
For the Thirring model, the non-observation of symmetry breaking is a robust result suggesting that the critical number of flavors required for dynamical symmetry breaking for the action (2) satisfies N c < 2. If we make the additional assumption that the UV fixed point of the Thirring model coincides with the IR fixed point of QED 3 , then this is compatible with the recent non-observation of a bilinear condensate in massless QED 3 with N ≥ 2 [28]. Beyond this, by contrast, the picture is not satisfactory. Despite the inapplicability of Elitzur's theorem, phase fluctuations due to the use of smeared sources have precluded fermion spectroscopy, and small fermion masses coupled with the absence of confinement have also prevented success in meson channels. Neither failure invalidates the DWF approach; the latter is simply a problem intrinsic to studying near-conformal physics on a finite volume, while the former might in future be tackled by a form of gaugefixing (see below). However, our inability to extract spectral information severely curtails insight into the failure of the axial Ward identity shown in Fig. 14, which potentially is a more profound problem, though again not necessarily fatal. What is disappointing, though, is that there is still no clear guide to the optimal lattice Thirring formulation.
Obvious future directions to explore include different formulations of the lattice Thirring model which may be less prone to the issues encountered here. Ref. [15] highlighted the role of a "hidden local symmetry" which is manifest once the Thirring action (13) is supplemented by a scalar Stückelberg field φ coupled to A µ . The HLS model has a gauge symmetry for which the Thirring model is the result of gauge-fixing to φ = 0.
Smoother gauge choices may enable better control over the fermion propagator. Formulating DWF with a true gauge symmetry will also restore transversity of the vacuum polarisation, making identification of the strong-coupling limit less ambiguous. Another route to finding and studying critical behaviour might come via implemention of the RHMC algorithm enabling N = 1 to be simulated. However, a more promising route, permitted by the control offered by the DWF formulation, may be to introduce a U(2N) and parity-invariant "Haldane" interaction term (ψγ 3 γ 5 ψ) 2 , motivated by the findings of the functional renormalisation group [4] which identifies a significant Haldane component in the fixed-point action corresponding to the Thirring model. There is still much to learn about fermions in 2+1d. andp µ and b(p) defined in (22). The hermitian operator D 0 D † 0 has zeromodes of the form ψ(s) = e ±αs : D 0 D † 0 ψ(s) = [b 2 (p) +p 2 − 2b cosh α(p) + 1]ψ(s), so the zero eigenvalue condition gives the definition of α in (22). The Green function of D 0 D † 0 is given by G 0 (s, s ′ ) = e −α|s−s ′ | 2b sinh α ≡ Be −α|s−s ′ | .
The solutions (55), (56) remain valid with m ↔ m 3 and ∆ ↔ ∆ 3 , but now Two features are apparent: first, there are no O(e −αLs ) contributions to ∆ 3 , so the first correction is O(e −2αLs ); second, the O(e −αLs ) contribution to A m 3 is now shifted by a phase e i π 2 with respect to the leading order piece. Both features mitigate finite-L s corrections to the calculation of ψ γ 3 ψ in the large-N GN model presented in Sec. 3.