Critical Behaviour in the Single Flavor Thirring Model in 2+1d

Results of a lattice field theory simulation of the single-flavor Thirring model in 2+1 spacetime dimensions are presented. The lattice model is formulated using domain wall fermions as a means to recover the correct U(2) symmetries of the continuum model in the limit where wall separation $L_s\to\infty$. Simulations on $12^3, 16^3\times L_s$, varying self-interaction strength $g^2$ and bare mass $m$ are performed with $L_s = 8, \ldots 48$, and the results for the bilinear condensate $\langle\bar\psi\psi\rangle$ fitted to a model equation of state assuming a U(2)$\to$U(1)$\otimes$U(1) symmetry-breaking phase transition at a critical $g_c^2$. First estimates for $g^{-2}a$ and critical exponents are presented, showing small but significant departures from mean-field values. The results confirm that a symmetry-breaking transition does exist and therefore the critical number of flavors for the Thirring model $N_c>1$. Results for both condensate and associated susceptibility are also obtained in the broken phase on $16^3\times48$, suggesting that here the $L_s\to\infty$ extrapolation is not yet under control. We also present results obtained with the associated 2+1$d$ truncated overlap operator DOL demonstrating exponential localisation, a necessary condition for the recovery of U(2) global symmetry, but that recovery of the Ginsparg-Wilson condition as $L_s\to\infty$ is extremely slow in the broken phase.


Introduction
The Thirring model is a quantum field theory of relativistic fermions interacting via a contact between conserved currents. In this paper we will examine the model in 2+1d spacetime dimensions, and therefore further specify the use of reducible (ie. fourcomponent) spinor fields, permitting the formulation of a parity-invariant mass term. The Lagrangian density reads where µ = 0, 1, 2 and a sum over i indexing N fermion species is implied. For a reducible representation of the Dirac algebra, there are two independent matrices γ 3 and γ 5 which anticommute with the kinetic term in (1), which results in a U(2N ) global symmetry generated by the set {1 1, γ 3 , γ 5 , iγ 3 γ 5 }. For m = 0 this breaks as U(2N ) →U(N )⊗U(N ). The symmetry can also break spontaneously through generation of a bilinear condensate ψ ψ = 0; while there are clear analogies with "chiral" symmetry breaking, this nomenclature should be eschewed in odd spacetime dimension. Spontaneous U(2N ) symmetry breaking is believed to be theoretically possible for sufficiently large self-coupling g 2 and sufficiently small N . Previous investigations of this question have used truncated Schwinger-Dyson equations [1,2,3], and Functional Renormalisation Group [4,5,6]. The prototype scenario has a symmetry breaking transition at g 2 c (N ), where g 2 c is a increasing function of N . A UV-stable renormalisation group fixed point can be defined as g 2 → g 2 c , where we identify a Quantum Critical Point (QCP) such that there exists an interacting continuum field theory solely specified by the field content, dimensionality and pattern of symmetry breaking. The critical flavor number N c required for symmetry breaking defined by g 2 c (N )| N =Nc → ∞, which in principle need not be integer, is an important property of the model. Since symmetry breaking is not accessible via expansion in any small parameter, the identification of N c for the Thirring model is an important and challenging problem in non-perturbative quantum field theory.
The model also provides a natural arena for the application of lattice field theory methods, and may well be the simplest fermion field theory requiring a computational solution. It turns out the choice of fermion discretisation has undue influence. Many early studies [7,8,9,10,11,12] used staggered fermions, which naturally support a symmetry breaking U( N 2 )⊗U( N 2 ) →U( N 2 ) [13]. The results of these studies are broadly consistent; the most wide-ranging of them, exploiting plausible assumptions about the g 2 → ∞ limit on a lattice, found N c = 6.6(1) [11]. The critical exponents found for N < N c appear rather sensitive to N , suggesting the existence of a rich family of distinct QCPs.
There are reasons to question whether this prediction for N c is correct; issues arising from a lattice perspective are reviewed in [14]. Here we present a non-rigorous plausibility argument for caution based on symmetry. The Thirring model in 2+1d may be studied analytically in the limit of large N , and the mass M V of a vector ff bound state interpolated by the current densityψγ µ ψ is predicted in this limit to be [15] M V m = 6π mg 2 ; lim Hence in the strong-coupling limit the Thirring model is a theory of conserved currents interacting via exchange of a massless spin-1 boson. Asymptotically the boson propagator switches from the canonical D µν (k) ∝ k −2 to ∝ k −1 ; this UV behaviour is exactly that predicted for the photon of QED 3 in the IR limit. Hence the Thirring QCP, if it exists, could well be identical with the IR fixed point of QED 3 , and the critical N c for both theories should then coincide. The parallels between the Thirring model and abelian gauge theory are more apparent still once an auxiliary vector field A µ is introduced, as reviewed in Sec. 2 below and more generally throughout the rest of the paper. Now, there is an old argument [16] for estimating N c in QED 3 , based on the conjecture and F is the thermodynamic free energy density. For asymptotically-free theories such as QED 3 f U V is related to the count of non-interacting constituents: Here 3 4 is the appropriate factor for Fermi-Dirac statistics in 2+1d, 4 is the number of spinor components per flavor and 1 counts the single physical polarisation state of the photon, which remains unconfined in QED 3 . For a phase with spontaneously-broken symmetry the number of weakly-interacting degrees of freedom is the Goldstone count plus the photon. Hence For QED 3 , and by extension the continuum Thirring model, the conjecture therefore predicts N c ≤ 3 2 , whereas for the symmetry breaking dictated by the staggered formulation the equivalent bound is the much less stringent N c ≤ 12. This disparity is a strong motivation for exploring lattice fermions with the correct global symmetry.
Recently this programme has been developed in two distinct directions. Lattice fermions employing the SLAC derivative, which is non-local but manifestly U(2N )symmetric, have been used to investigate the model in [17,18]. Since the Thirring model is not a gauge theory, the long-standing objections [19] to the SLAC approach associated with lack of localisation of the fermion-gauge vertex do not apply. No evidence has been found for spontaneous symmetry breaking for any integer value of N (ie. any unitary local version of the model); an estimate N c = 0.80(4) < 1 was reported in [18]. Meanwhile, the Thirring model formulated with domain wall fermions (DWF) has been investigated by us in [20,14]. DWF employ a local lattice derivative at the cost of introducing a fictitious extra dimension. In this case the U(2N ) symmetry is not manifest but hopefully recovered in a controlled way in the limit that the separation L s between domain walls located at either end of this "third" dimension is made large. An aspect worth highlighting is that in the most promising approach the fermions interact with the auxiliary field A µ throughout the bulk, ie for 0 < x 3 < L s , very similar to the way gauge theories are formulated with DWF [21].
We have found that the L s → ∞ limit becomes particularly challenging precisely in the strong-coupling regime where symmetry breaking is anticipated. In particular, Ref. [14] presented results for N = 0, 1, 2 on 12 3 × L s (for N > 0) with L s ranging from 8 to 40. The results for N = 0 and N = 2 are fairly clear-cut: symmetry-breaking is observed in the former case and not in the latter. For N = 1 there are qualitative indications of a change in the system's behaviour at the strongest coupling explored (g −2 a = 0.3, where the lattice spacing a defines the scale). In particular the bilinear condensate ψ ψ(m) displays significant L s -dependence in this regime, and the resulting signal estimated in the L s → ∞ limit is significantly greater than that from weaker couplings; there is also a marked departure from linear m-dependence. The conclusion reached in [14] is that 0 < N c < 2 with "strong evidence" for N c > 1.
Settling the issue is hard for a couple of reasons. One is that the RHMC simulation algorithm required for N = 1 is numerically more demanding [14]. Another is that in a finite volume the bilinear condensate vanishes identically for m = 0 so that neither order parameter nor associated susceptibility are directly accessible in the U(2N )-symmetric limit. In the current paper we apply improved code and substantially enhanced computing resources to both issues, by exploring the model on both 12 3 × L s and 16 3 × L s , with L s as large as 48, with much finer resolution along the coupling axis, particularly in the suspected critical region. The resulting order parameter data ψ ψ(g 2 , m) are obtained with sufficient statistical precision to permit fits to a renormalisation-group inspired equation of state (EoS) applicable away from m = 0, yielding estimates for both critical coupling g 2 and accompanying critical exponents β m and δ defined in (15) below. A similar strategy has been applied successfully in staggered fermion studies [7,9,11]. As the methodology implies, the main results of the study are confirmation that N c > 1 and a first tentative characterisation of the critical properties of the QCP.
The remainder of the paper is organised as follows. Sec. 2 sets out the lattice formulation of the model using DWF and the auxiliary field method, and reviews the simulation algorithm and principal observables. Results are presented in Sec. 3 in three subsections. Sec. 3.1 presents a systematic study of varying β ≡ g −2 a, m, L s and space-time volume, in the regime 0.3 ≤ β ≤ 0.52. Results for ψ ψ are first extrapolated to L s → ∞, and then fitted to the empirical EoS (16) below. The picture that emerges appears under control; the L s -extrapolation and EoS-fitting procedures almost, but not quite, commute, yielding fairly stable estimates for β c ≈ 0.28 and estimates for the exponents β m , δ distinct from those of mean field theory. Finite volume effects appear remarkably small. However, all the data used in this analysis lie in the symmetric phase. To address this Sec. 3.2 presents a complementary study at fixed L s = 48, but this time exploring couplings as strong as β = 0.23. Both ψ ψ and its associated susceptibility χ are studied. Here some issues emerge; the EoS fit is not so successful at describing data in the strong-coupling phase, and χ , while having a peak as expected near criticality, exhibits a non-standard scaling as m is varied. We also present results for a residual δ h introduced in [22] to quantify recovery of U(2N ) symmetry, and the bose action density (2g 2 ) −2 A 2 µ . In Sec. 3.3 we present results for the locality, and recovery of the Ginsparg-Wilson relation, of the truncated overlap operator (corresponding to finite L s ) in the critical region. Both are key ingredients in demonstrating the existence of a local U(2N )-symmetric field theory at the QCP. Our conclusions are discussed in Sec. 4.

Formulation and Numerical Simulation
We use an auxiliary field formulation to represent the Thirring interaction, which recasts the fermion action as a bilinear form while preserving the global symmetries. In continuum notation the Lagrangian density reads On a lattice the vector auxiliary field A µ is naturally formulated on a link. Preserving fermion global symmetries while transcribing to a lattice is of course a long-standing issue in lattice field theory. Our approach uses domain wall fermions (DWF), in which a fictitious third dimension is introduced so the fermions Ψ(x, s) are defined in 2+1+1d: We use the Möbius formulation, implying that D W [A = 0], D 3 are free Wilson derivative operators, with D 3s,s having open boundary conditions implemented on the hopping terms, viz. δ s∓1,s (1 − δ s ,1/Ls ), at domain walls located at s = 1, L s . The auxiliary field A µ (x) is 3-static, taking the same value on every spacetime slice along the 3rd direction.
In previous work we have referred to this as the bulk version of the model [14]. Our model employs a non-compact formulation of the interaction in which each hopping term in D W carries a non-unitary link factor (1 ± iA µ (x)); in this way integration over A generates solely 4-fermi terms, and not higher point contact interactions. Further details are set out in [14].
In the large-L s limit the free kinetic operator approaches an overlap operator D OL defined in 2+1d and satisfying a generalisation of the Ginsparg-Wilson relations [23,24] For weakly-interacting fields, the RHS of the first two of these relations are formally O(a) and thus vanish in the continuum limit, recovering the desired U(2) global symmetry of the continuum model. This only holds for the bulk formulation of the Thirring model with DWF, and not the alternative "surface" formulation discussed in [20,14]. For strongly-interacting fields the recovery of the GW relations (8) and the locality of the corresponding D OL will be discussed further in Sec. 3.3 below. The bridge between 2+1+1d and the target model rests on the identification of physical fermion fields in 2+1d localised entirely on the domain walls, which we regard as a working assumption: with projectors P ± ≡ 1 2 (1 ± γ 3 ). The mass term mS m in (7) needs some discussion. A U(2)-symmetric theory has three physically equivalent parity-invariant mass terms m hψ ψ, im 3ψ γ 3 ψ and im 5ψ γ 5 ψ. For DWF fermions with L s finite the choice imψγ 3 ψ with ψ,ψ specified in (9) gives the best approach to L s → ∞ [22,24,14], and we use this definition throughout this paper. Approach to the U(2) symmetric limit will be monitored in Sec. 3.2 below via the residual δ h ψ ψ − ψ iγ 3 ψ [22]. Specifying the bilinear component of the action (7) by M, then the action simulated using bosonic pseudofermion fields Φ, Φ † is Assuming detM is real, then the resulting functional weight det[M m 3 =m M −1 m h =1 ] tends to detD OL (m) in the limit L s → ∞ [24]. An RHMC algorithm is needed to handle the fractional powers in (10), as described in [14]. Some tests of the accuracy of the rational approximation needed to implement fractional powers in the parameter regime of interest are presented in Sec. 3.2. In the meantime the code has been modified in two main aspects: multiprocess parallelisation and a simplified, more stable solver. The parallelisation makes use of the MPI paradigm, with a 3D domain decomposition across the spatial and temporal directions, leaving the Domain Wall dimension uncut to allow the compiler -as in the original version of the code -to perform automatic vectorisation. The alternative solver implementation is a conjugate gradient-based multi-shift solver, having the advantages of requiring 33% less memory compared to the original QMRbased one, and being stable in single precision. The possibility of running the solver in single precision during the molecular dynamics part of the algorithm leads to another 50% save in memory which, as the solver is severely memory-bound, translates into a direct performance increase. We have made the simulation code publicly available [25].
Finally, the observables studied in the bulk of the paper are simply the bilinear condensate and its associated susceptibility The notation is slightly loose; the bilinear actually used in computations should be understood to be iψγ 3 ψ. The susceptibility χ defined in (12) is often referred to as the disconnected component, and is expected to manifest any singular behaviour expected at a continuous phase transition. The full susceptibility corresponding to ∂ 2 ln Z/∂m 2 contains an extra connected component not included in the variance of the bilinear order parameter calculated here. In our work the expectations (11,12) are calculated with unbiassed estimators using 10 independent Gaussian noise vectors per configuration. First let's discuss the L s → ∞ extrapolation; empirically the convergence to the large-L s limit is exponential [14]: Sample fits are shown in Fig. 1. The extracted decay constant ∆(β, m) is shown in Fig. 2. While errors are appreciable, particularly at weaker couplings where the signal is small, it can be seen that for the weaker couplings 0.36 ≤ β ≤ 0.44 ∆ ≈ 0.03 and is approximately m-independent, whereas for the stronger couplings 0.3 ≤ β ≤ 0.34, ∆(m) is roughly linear, with an intercept ∆(m = 0) ≈ 0.007. This suggests that L s ∼ O(10 2 ) is needed to completely control the extrapolation in this regime.
Next we turn to analysis of the critical point, assumed present at (β, m) = (β c , 0). Since data are collected in the presence of a U(2)-symmetry breaking mass m = 0, we follow an indirect route, motivated by the renormalisation group [7]. At fixed spacetime volume, assume a scaling form where t ≡ β − β c . In the limits m = 0, t = 0 we immediately recover whereupon we identify the conventional critical exponents β m and δ familiar from the ferromagnetic transition. For small t we may Taylor-expand the scaling function F to yield the equation of state (EoS):  Results of a least-squares fit of (16) to data from 16 3 × 48, 0.3 ≤ β ≤ 0.52, and ma = 0.005,0.01,0.02,. . . , 0.05 are shown in Fig. 3, together with the curve resulting from the same fit parameters as m → 0. The data support a symmetry-breaking continuous phase transition at β c = 0.2537(2) (the fit of Fig. 7 below with significantly smaller χ 2 /dof results from excluding the data with β = 0.3; here we retain them for the sake of uniformity in the subsequent analysis). The fitted exponents β m = 0.42(1), δ = 3.41(5) differ significantly from their predicted values in mean field theory, namely β m = 1 2 , δ = 3.
However, in order to identify this transition with the breaking of U(2) symmetry a minimum requirement is stability under the extrapolation L s → ∞. The purist's approach, presented first, is to fit (16) to data which has been first extrapolated using (13). Such a fit is shown in Fig. 4. The first thing to note is the far larger errors on the datapoints (and correspondingly smaller χ 2 /dof) -the extrapolation is particularly difficult to control at weak coupling where the signal is small. The fitted EoS still supports a continuous phase transition, though with modified exponents β m = 0.31(2), δ = 4.3 (2). The critical point is also shifted to weaker coupling: β c = 0.279(1). In the immediate vicinity of the critical point the fitted curves do not well describe the small-mass data ma = 0.005, 0.01; however, due to the large errorbars exclusion of these masses does not significantly alter the fit.
A more pragmatic approach is to perform fits at fixed L s and then study the be- haviour of the fit parameters as L s → ∞. This is less well-motivated theoretically, but has the practical advantage that the data passed to the least-squares fitting procedure is of much higher statistical quality. Results for the critical coupling and the exponents are shown in Figs. 5, 6 respectively. In Fig. 5 it is notable that β c (L s ) from 12 3 and 16 3 volumes are compatible; the evolution with L s is smooth but falls significantly short of the fit of the extrapolated data even by L s = 48, so it is not clear whether extrapolation and fitting commute. In Fig. 6 by contrast the fitted β m (L s ) and δ(L s ), while agreeing at large L s , approach this limit from opposite directions on 12 3 and 16 3 ; moreover while on 12 3 the data from the largest available L s are compatible with the values extracted from the extrapolated data, as already remarked there is a significant disparity on 16 3 .

Towards Stronger Couplings on 16 3 , L s = 48
The results of the previous subsection support the claim made in [14] that the N = 1 model has a continuous phase transition associated with the onset of a bilinear fermion condensate for β ≈ 0.3. There is no sign of any significant finite volume effect. Moreover, for the largest L s examined, the transition is reasonably well-modelled by an EoS with critical exponents both distinct from mean-field values, and consistent with a QCP defining a previously unknown strongly coupled quantum field theory of fermions. We might be concerned, however, that all data analysed in Sec. 3.1 lie on the weak-coupling side of the putative transition, and that there are hints from Fig. 3 and particularly   To investigate further we now turn to data generated using approximately 5000 RHMC trajectories at fixed L s = 48 on spacetime volume 16 3 , but extending to β-values on the strong-coupling side of the transition. Additionally, at the lightest mass ma = 0.005 the β-axis is sampled with a finer resolution ∆β = 0.01 in the strong coupling region. Fig. 7 shows the ψ ψ(β, m) dataset, and a fit to (16) based on 0.32 ≤ β ≤ 0.52. The fit is compatible with that shown in Fig. 3, and as advertised is of slightly higher quality as a result of excluding β = 0.30. It is immediately apparent, however, that it fails to model data on the strong coupling side of the transition; here ψ ψ falls below the model, the effect is more pronounced with decreasing m, until by ma = 0.005 the curve becomes flat for β < ∼ 0.25, as shown in Fig. 8. The flattening of the condensate at strong coupling is not a new story; indeed, simulations with staggered fermions actually exhibit a maximum before dropping as β → 0 + . A possible origin of this behaviour was suggested in [7]: as a result of the linear coupling between A µ and the fermion current, the leading order large-N correction to the auxiliary propagator is not transverse in a lattice regularisation, leading to breakdown of positivity as β → 0. In the large-N approach the effect is mitigated by an additive renormalisation of g −2 , so that the strong coupling limit of the model is now taken as β → β * > 0. In [11]   consequent prediction of N c . A decrease of ψ ψ β→0 has also been reported in simulations with SLAC fermions [17,18] and with DWF in a variant "surface" formulation of the model [20], suggesting that strong coupling lattice artifacts are a generic feature of the Thirring model, and may have a more general origin than that suggested by the large-N approach. Be that as it may it will clearly be important to establish a clear separation between β * and any β c associated with a Thirring model QCP. From Fig. 8 we might estimate β * ≈ 0.25, uncomfortably close, with current resolution, to the β c estimates of Sec. 3.1. At this point it is appropriate to discuss a technical aside. In the RHMC algorithm described in [14], it is necessary to calculate fractional powers of the fermion kernel A = M † M. In practice this is performed using a rational approximation where the coefficients α i , β i may be calculated using the Remez algorithm implementation described in [26]. They are chosen so that over a spectral range (λ d , 50.0), |r p (x) − x p | < 10 −6 for matrices needed during trajectory guidance and < 10 −13 for those needed in the Monte Carlo acceptance step. For all work to date we have used λ d = 10 −4 corresponding to the smallest value of (ma) 2 explored, which translates to partial fraction numbers N pf = 12 (guidance) and N pf = 25 (acceptance); however one might question whether this is sufficiently accurate for studies with ma = 0.005. Accordingly we have performed "enhanced" simulations at three β-values with Remez coefficients generated with λ d = 10 −5 , corresponding to N pf = 14 (guidance), and N pf = 29 (acceptance). As shown in Fig. 8, fortunately there appears to be no significant difference with data calculated using the previous λ d = 10 −4 . Next we present data for the susceptibility χ defined in (12), for the whole dataset in Fig. 9 and for the lightest ma = 0.005 in Fig. 10. As might be anticipated, statistical errors in χ are considerably larger than those for the condensate, and accordingly we choose not to attempt an L s → ∞ extrapolation. However, again, the agreement between results obtained using the default and enhanced rational approximations seen in Fig. 10 is reassuring. For each value of m χ (β) is non-monotonic, the peak shifting to stronger coupling as m decreases in accord with expectations for a second derivative of the free energy at a critical point; this is corroborated by the model prediction obtained by differentiation of (16) with respect to m, and plotted using the fitted parameters in Fig. 11. Fig. 10 suggests that the location and even the sharpness of the peak at criticality is roughly as expected, once an empirical rescaling is applied. However, there are features of Fig. 9 which are clearly problematic; the m-ordering of the data is opposite to model expectations, with χ increasing with m over the whole β-range studied, and the convergence of χ curves with different m as β grows large seen in Fig. 11 is not observed. We postpone further discussion of these issues to Sec. 4.
Next we discuss the approach to recovery of U(2) symmetry expected as L s → ∞. Ref. [22] introduced a residual δ h defined in terms of the 2+1+1 dimensional fields as follows: In [22] 2δ h was found to furnish a lower bound for the difference between ψ ψ and ψ iγ 3 ψ , and to vanish ∝ e −cLs for quenched QED 3 . In the Thirring model, δ h (L s ) for various couplings was presented in Fig. 12 of [14]. While in all cases δ h still decreases with L s , it grows larger as coupling increases, and by β = 0.3 its decay constant c even develops a dependence on m. Fig. 12 taken at fixed L s = 48 confirms that m-dependence of δ h does indeed set in for β < ∼ 0.4, and that δ h continues to grow as β decreases, suggesting that recovery of U(2) symmetry will be an ever-increasing challenge in the symmetry broken phase as m → 0.
Finally, Fig. 13 shows results for the bosonic auxiliary action density (2g 2 ) −1 A 2 µ . As discussed in [14], for DWF with the bulk formulation of the Thirring model there is no simple interpretation in terms of a local four-fermion condensate available; rather we regard it as an extra observable sensitive to light fermion dynamics. Its behaviour is non-monotonic, with a minimum at β 0.46 before rising to approach and then exceed the free-field value 3 2 at β 0.24. The notable feature of Fig. 13 is the fermion massdependence; broadly speaking the departure from the free-field result increases with decreasing m (although the m-ordering of the data is somewhat noisy), the effect being most pronounced for 0.3 < ∼ β < ∼ 0.4 immediately above the suspected critical region.

Properties of the Associated Overlap Operator
The equivalence of DWF [27,21] and the (truncated) overlap operator [28] is well established in 3+1d, eg. [29]. This equivalence is further shown in 2+1d [24] for both the regular mass term mψψ and the linearly independent twisted mass terms imψγ 3,5 ψ introduced above. As such, locality of the domain wall operator in the target dimensionality can be demonstrated by showing the locality of the overlap operator.
We use the Shamir and Wilson formulations of the overlap operator with twisted mass −imψγ 3 ψ given by where the Wilson and Shamir kernels, defined via V = γ 3 sgn(H), are and D W ≡ D W (−M ) is the massive Wilson Dirac operator, M fixed to 1. Note, however, that the link factors multiplying the difference operators in D W are the non-unitary [1 ± iA µ ] rather than the unitary e ±iAµ characteristic of a gauge theory. The key relation Hγ 3 = (γ 3 H) † is preserved. Our formulation of DWF in the L s → ∞ limit is expected to recover (19) with the Shamir kernel [24].
The standard mass formulation of the overlap operator is D I OL (m) = 1+m 2 + 1−m 2 V S/W . The signum function typically is approximated with a rational function resulting in a truncated overlap operator. The hyperbolic tangent (polar) approximation is the most commonly used and may be expressed as [30] sgn(x) ≈ tanh(ntanh −1 x) = xn for n even. For the Shamir formulation it is much more efficient to use a formulation exploiting an extra dimension, and the truncated overlap operator can be reconstructed directly from a domain wall formulation. The standard mass and alternative mass domain wall operators may be expressed as D Then with we have the following relation [24,29], where the precise form of the terms is unimportant for our purposes.
So the overlap operator (19), with the Shamir kernel (20), truncated with the polar approximation (21), evaluated with given n, is identical to the top left entry of matrix K (25) using domain wall extent L s = n. There is a similar relation for the Wilson kernel, with a different domain wall formulation. We also have K I = C † (D I DW (1)) −1 D I DW (m)C.
Invariant transformations, corresponding to the Ginsparg-Wilson (GW) relations (8) are given by [22] Ψ This relation is exact for the 2+1d overlap operator D = D OL , and is reproduced by DWF in 2+1+1d in the L s → ∞ limit. In order to recover the U (2) symmetry in the continuum limit a → 0, we must have the GW terms aDγ 3,5 D and equivalently the transform terms aD 2 vanishing in the same limit. A sufficient condition for this to be the case is the Dirac operator being exponentially local. Evidence for this has been given for the overlap operator in 3+1d [31], and it behooves us to investigate the 2+1d case. Since the overlap operator is a dense matrix and manifestly non-local, demonstration of exponential locality is important, especially around critical regions.
To see that exponential locality is necessary to ensure recovery of the continuum U (2) symmetry, note that so that recovery requires [aDΨ] a→0 = 0 (28) < ∞ which is true if D is exponentially local, and hence exponential locality allows recovery of U (2) symmetry. In order to illustrate the locality of the overlap operator, we follow [31]. Let where the point source η(x) = δ x,y δ α,1 where y is an arbitrary location and α a spinor index. Then we calculate the decay as where the "Manhattan taxi distance", ||x − y|| 1 = µ |x µ − y µ | is just the L 1 norm. The locality of D OL in the critical region is illustrated for the twisted imψγ 3 ψ mass term with ma = 0.005 in Figs. 14,15. Within the limitations imposed by a 16 3 volume, Fig. 14 is consistent with exponential falloff for ra > ∼ 10. There is a mild β-dependence; as expected the falloff slows down as the coupling gets stronger. Also shown for comparison are data obtained with unitary link fields e ±iAµ , where exponential localisation is more manifest. Convergence to a meaningful value of the decay rate on this small volume is seen to be difficult in Fig. 15 showing the decrement from one r-value to the next, by plotting f (r)/f (r − 1). As the coupling strength gets closer to criticality there is only   We also examine the truncation/finite-L s error of the overlap/DWF operator via the GW term, as a means to assess recovery of U(2) symmetry. In the n → ∞ (L s → ∞) limits the GW error δ GW , given as with φ a complex field chosen at random at each lattice site such that each component is distributed uniformly in (−1, 1), should be exactly zero for zero mass. We use a single auxiliary field configuration to plot δ GW in Fig. 16. Although the boson field is generated with a non-zero mass, the GW error is measured with D(m = 0). The GW error vanishes only very slowly, and in fact δ GW is numerically very similar to δ h defined in (18) (plotted as a function of L s in Fig. 12 of [14]). It is interesting that convergence with the Wilson kernel is slightly faster as compared to the Shamir formulation, which is surprising since with non-unitary links the Wilson kernel is not bounded which should prejudice convergence [30,24]. Again, for comparison results obtained with unitary link fields are included in the inset; here exponential falloff of the error is much sharper. These results suggest that the very large values of L s needed for U(2N ) recovery in the critical region have their origin in the non-unitary nature of the link fields in this formalism.

Discussion
The main conclusion of this study is that we can state with much more confidence that there is a phase transition associated with bilinear condensation in the Thirring model defined with N = 1 domain wall fermions, and hence that N c > 1 as originally suggested in [14]. We have also made the first steps towards characterising the critical properties. The enhanced dataset we have generated demonstrates both the importance and the stability of the L s → ∞ extrapolation, enabling fits to an RG-inspired equation of state with β c 0.279(1) which work extremely well for β > ∼ β c . The fitted critical exponents are significantly different from their mean-field values, although both Fig. 2 and Figs. 5,6 suggest simulations with larger L s , perhaps even L s ∼ O(10 2 ), will be needed in order to fully control the predictions. Further encouraging signs are a susceptibility peak of the expected shape seen in Fig. 10 and the m-dependence of the bosonic action in the subcritical region seen in Fig. 13.
The result N c > ∼ 1 is in clear contast with predictions obtained with staggered fermions, for the reasons reviewed in the Introduction, but more interestingly also with N c = 0.80(4) obtained with SLAC fermions [18]. We speculate that the two lattice approaches describe different continuum theories, and that the bulk DWF formulation followed here more closely conforms to a picture of the strong dynamics in which the auxiliary boson A µ resembles a gauge field. Ultimately, physics at a QCP is specified not by a Lagrangian density, either continuum-like or regularised, but by more primitive considerations such as dimensionality, field content, and of course the pattern of global symmetry breaking. Interestingly, a recent study using the Conformal Bootstrap predicts 1 < N c < 2 for QED 3 [32].
For the first time in Sec. 3.3 we have studied properties of the 2 + 1d overlap operator D OL associated with DWF in the vicinity of the critical point. Fig. 14 suggests there exists a limit where aD OL → 0, a necessary condition for the existence of a continuum limit with conventional U(2N ) symmetry, although confirmation will require studies on larger spacetime volumes. However, symmetry recovery also requires the restoration of the Ginsparg-Wilson relations (8). Both the direct estimates shown in Fig. 16 and the indirect measure via the residual δ h shown in Fig. 12 suggest this is at best restored only very slowly in the critical region. Note that δ h → 0 is not guaranteed even once O(e −∆Ls ) corrections are applied to the order parameter as described in Sec. 3.1. It is now becoming apparent that this behaviour has its origins in the use of non-unitary link fields in the Wilson operator D W ; recall this originates in the desire to have only four-point fermion interactions left once the auxiliary is integrated over. Fig. 16 is also a motivation to explore measurement using D OL defined with the alternative Wilson kernel, and/or with expansion order n greater than the L s used in ensemble generation with DWF.
It's also important to review areas where the simulation falls short of the QCP ideal. As outlined in Sec. 3.2, at fixed L s = 48 EoS fits fail to describe the data with β < β c , and the susceptibility plot Fig. 9 shows an inverted m-hierarchy when compared with the expectation of Fig. 11, and moreover requires a seemingly arbitrary rescaling even to fit a single m value. The behaviour at weaker couplings can be somewhat accommodated by replacing the LHS of the EoS (16) by a factor m(m/m 0 ) α . The inset of Fig. 11 shows the susceptibility curves thus obtained with m 0 a = 0.005, α = 0.3. This fails to fix the hierarchy in the critical region, however; it seems safer to conclude that χ defined in (12) is not the second derivative of the free energy of a 2+1d theory. A modification of the physical field prescription (9) may be needed, to allow for the possibility that the relevant modes close to the walls "leak" somewhat into the bulk as ma → 0. A related puzzle, again a failure to reconcile the observed behaviour of order parameter and susceptibility data with possibly the same cause, is the breakdown of the axial Ward Identity noted in [20].
The flattening of the order parameter at small β seen in Figs. 7,8 remains a worrying aspect. Possibly the L s → ∞ extrapolation is not yet under control in the broken region -this scenario is supported by the results of pilot simulations on 16 3 × 64 with β = 0.3, ma = 0.005, plotted in Figs. 8,10 (the L s = 64 point lying on the fitted curve in Fig. 8 should be regarded at this stage as a coincidence). Another possibility is that the scaling window where (16) is applicable is simply very narrow on the horizontal scale used in these figures. As things stand, however, we have not yet demonstrated a clear separation between β c and a putative β * where lattice artifacts dominate, and hence do not yet have an understanding of the broken phase comparable with that for β > β c . The slow convergence to the U(2N ) limit discussed above may prove a formidable obstacle; there is much work still to be done.