Domain Wall Fermions for Planar Physics

In 2+1 dimensions, Dirac fermions in reducible, i.e. four-component representations of the spinor algebra form the basis of many interesting model field theories and effective descriptions of condensed matter phenomena. This paper explores lattice formulations which preserve the global U(2N ) symmetry present in the massless limit, and its breakdown to U(N)xU(N) implemented by three independent and parity-invariant fermion mass terms. I set out generalisations of the Ginsparg-Wilson relation, leading to a formulation of an overlap operator, and explore the remnants of the global symmetries which depart from the continuum form by terms of order of the lattice spacing. I also define a domain wall formulation in 2+1+1d, and present numerical evidence, in the form of bilinear condensate and meson correlator calculations in quenched non-compact QED using reformulations of all three mass terms, to show that U(2N) symmetry is recovered in the limit that the domain-wall separation tends to infinity. The possibility that overlap and domain wall formulations of reducible fermions may coincide only in the continuum limit is discussed.


Introduction
Relativistic fermions moving in two spatial dimensions are at the heart of many interesting issues in theoretical physics, and are increasingly important in effective descriptions of phenomena in condensed matter systems. Let us immediately draw a distinction between so-called irreducible formulations in which the fermion fields ψ,ψ are described by two-component spinors, and reducible models which invoke four-component spinors. In the former case, there is no mass term invariant under a discrete parity inversion; gauge theories of irreducible fermions generically manifest a parity anomaly [1,2,3], leading to an induced Chern-Simons term endowing the gauge degrees of freedom with mass. Such theories describe excitations with fractional statistics, and form the basis for effective descriptions of the fractional quantum Hall effect [4]. The main focus of this paper, by contrast, will be reducible theories, which like their 4d counterparts admit parityinvariant fermion masses, but whose action is invariant under an enlarged group of global symmetries generated by both γ 5 and the "unused" γ 3 . The prototype model involving reducible fermions is QED 3 , an asymptotically-free theory with potential non-trivial critical dynamics in the infra-red, and still investigated as a toy model of walking technicolor [5]. Reducible fermions figure in effective descriptions of the spin-liquid phase of quantum antiferromagnets [6,7], the pseudogap phase of cuprate superconductors [8,9], and graphene [10,11,12].
Once interactions are introduced, reducible 3d fermions may exhibit spontaneous mass generation via formation of a bilinear condensate ψ ψ = 0, entirely analogous to chiral symmetry breaking in QCD. This is manifested in variants of the Gross-Neveu (GN) model [13], in which the basic interaction is a four-fermi contact between scalar or pseudoscalar densities; and also in the Thirring model [14] whose interaction is a contact between two conserved currents, and which shares the global symmetries of QED 3 . For GN models the condensate may be studied systematically using an expansion in 1/N f where N f is the number of fermion flavors [15], whereas condensate formation in the Thirring model is inherently non-perturbative. In both cases the critical coupling at which the condensate forms is associated with a UV-stable fixed point of the renormalisation group, although the nature of the fixed-point theories remains an open problem [16].
Because field fluctuations near a fixed point can be large, it is natural to explore critical dynamics in 3d using lattice field theory techniques. We immediately confront the issue of how best to formulate the fermion fields. The two traditional approaches which mitigate the species doubling present in naive discretisations are Wilson fermions, which eliminate unwanted species by explicitly breaking important global symmetries, and staggered fermions, which preserve a subgroup of the symmetries by effectively spreading the spinor degrees of freedom over several lattice sites (for an excellent introduction to lattice fermions see chapters 5,7 and 10 of [17]). In a seminal work [18] Coste and Lüscher argued that the Wilson formulation is the more natural for irreducible descriptions, and can be shown to recover the correct parity anomaly (a recent calculation of the mass-dependence of the anomaly in various background fields using Wilson fermions is given in [19]), while the staggered approach naturally leads to a parity-invariant mass term; indeed the relation between staggered fields and the N f = 2 reducible spinor flavors expected in the continuum limit was given by Burden and Burkitt in [20]. Since the Chern-Simons action is imaginary in Euclidean metric, numerical simulations of 3d fermions have to date been restricted almost exclusively to this second case.
To date there have been several papers studying critical behaviour in both GN [15,21,22] and Thirring [23,24,25] models using numerical simulations of staggered fermions. In general the results support the theoretical prejudice that GN models exhibit a critical point at a coupling g 2 ∼ O(a) (where a is lattice spacing) for all N f , with small corrections of O(1/N f ) to both g 2 c and the "mean field" critical exponents predicted in the large-N f limit [15]. By contrast, the Thirring model exhibits gap formation only for N f ≤ N f c , with both exponents and g 2 c strongly dependent on N f [25]. In the continum, the two models are supposedly distinct. However, recent results obtained with a fermion bag algorithm permitting simulation in the massless limit [26,27] suggest that in fact the lattice models defined using staggered fermions may actually lie in the same universality class for the minimal case N f = 2. This somewhat unexpected result is motivation to question the applicability of staggered fermions to problems involving critical fields. For N staggered flavors the global symmetry is U(N)⊗U(N), with N f = 2N, broken to U(N) by the generation of a fermion mass. In the long wavelength limit a∂ → 0 this enlarges to the required U(4N) →U(2N)⊗U(2N) [20], but there is less reason than usual to trust this restoration once fluctuations on all length scales are important. Other situations where it may be important to reproduce the global symmetry pattern correctly include the infrared behaviour of QED 3 , where it has been argued that dynamical mass generation depends on the correct counting of massless degrees of freedom, which include Goldstone modes in a symmetry-broken phase [28], and the role of "half-instantons" in the thermal response of the gapped phase of graphene, which may result in a much reduced value for the critical temperature [29].
This paper explores the application of formulations originally developed to optimise the reproduction of global symmetries in lattice QCD, namely Ginsparg-Wilson (GW) fermions [30] and, principally, domain wall fermions [31,32], to reducible fermion models in 2+1d. After reviewing the relevant symmetries and identifying three distinct but physically equivalent formulations of the mass term in the next section, in Sec. 3 we generalise the GW relation to fermions in 2+1d and identify remnant quasi-global symmetries, which recover the desired U(2N f ) form only in the continuum limit a → 0. A realisation of the GW symmetries by an overlap operator [33] is given. In Sec. 4 we define a domain wall fermion operator in 2+1+1d which permits the definition of fermi fields localised on domain walls at either end of the newly introduced 3 direction which purport to satisfy the U(2N f ) symmetry in the limit that the wall separation L s → ∞. An important component of the argument is the reformulation of the three distinct mass terms given in Sec. 2. Sec. 5 presents results from numerical investigations of the N f = 1 domain wall operator in the context of quenched non-compact QED 3 , which permits the use of either weak, strong, or intermediate coupling. While there is no attempt to explore either continuum or thermodynamic limits, we calculate both bilinear condensates (Sec. 5.1) and meson propagators (Sec. 5.2) using each of the three alternative mass terms, and show that in almost all cases as L s → ∞ the results are in accord with a scenario in which U(2) symmetry is broken to U(1)⊗U(1). Interestingly, the most rapid convergence to the U(2)-symmetric limit is obtained for the case of a "twisted" mass term imψγ 3 ψ. For intermediate coupling the results for the condensate ψ ψ are compatible in the massless limit with old results obtained with staggered fermions [35]. Finally in Sec. 6 we present a summary of the findings and an outlook for future investigations. We also discuss the intriguing possibility that for reducible theories of fermions in 2+1d the overlap and domain wall approaches may not coincide except in the continuum limit.

Relativistic Fermions in 2+1d
I begin by reviewing the continuum formulation of a gauge theory with fermion fields Ψ,Ψ in a reducible representation of the spinor algebra, based on 4 × 4 Euclidean Dirac matrices γ µ with {γ µ , γ ν } = 2δ µν , µ, ν = 0, 1, 2, and having a parity-invariant mass. The weakly-interacting long-wavelength limit of staggered lattice fermions naturally reproduces this formulation with N f = 2 flavors [20] -in what follows flavor indices are suppressed. The action can be written (for convenience, the necessary d 3 x is omitted in all action definitions) where the covariant derivative operator D can be expanded as This has global symmetries where γ 3 and γ 5 are two additional traceless, hermitian, and linearly independent 4×4 matrices which anticommute with all the γ µ (see (8,9) below), and as usual in Euclidean matric γ 5 ≡ γ 0 γ 1 γ 2 γ 3 . For fermion mass m = 0 there are two additional symmetries These four rotations generate a global U(2) invariance, which generalises to U(2N f ) for several flavors. The mass term explicitly breaks the symmetry from U(2N f ) → U(N f )⊗U(N f ). It will prove interesting to explore different forms of the mass term, which are simply accessed by changing integration variables in the path integral. Since there is no axial anomaly in 2 + 1d, this procedure is straightforward in the continuum and the resulting action describes identical physics. If, however, the representations of the Dirac matrices are tied to the particular form of the underlying lattice, as is the case for staggered fermions or graphene, then due to discretisation effects the mass terms are not equivalent and correspond to distinct patterns of symmetry breaking (see the discussion following eqn. (17) for an example). Let's recast the continuum action (1) in terms of two twocomponent spinors u and d: whereD = −D † = σ 1 D 0 + σ 2 D 1 + σ 3 D 2 and the σ i are Pauli matrices. The link with (1) requires the identification implying We now define an important discrete symmetry, parity, here specified for convenience in terms of reversal of all three spacetime axes x µ → −x µ (in general parity must invert an odd number of axes, since flipping an even number is equivalent to a rotation: the Euclidean parity operation which flips just one axis is formally equivalent to the timereversal operation frequently discussed in condensed matter physics). In fact it can be realised in two ways: i.e. Ψ → iγ 5 Ψ;Ψ → −iΨγ 5 (11) This should be no surprise, since both γ 3 and γ 5 behave identically with respect to the γ µ appearing in (1). In either case the parity operation effectively exchanges the u and d fields, absorbing the sign change ofD under x → −x, but keeping the mass term invariant. Now consider a change of basis in which the action (7) reads The parity transformation leaving (13) invariant is equivalent to (10). In terms of fourcomponent spinors the action (13) can be written Although the mass term imΨγ 3 Ψ is still parity invariant, it is now an antihermitian term in the Lagrangian, in contrast to the hermitian mass of (1). Superficially it resembles the twisted mass sometimes used to formulate lattice QCD in 4d [36], though the absence of an anomaly in 2 + 1d permits the term to be flavor singlet. We can also consider a different change of basis yielding with this time a parity transformation (11). In terms of four-component spinors (16) can be written S =Ψ(D − imγ 5 )Ψ; (17) again, the mass term is parity invariant and antihermitian. While the three actions (1,14,17) must be equivalent in the continuum, when derived from a lattice system such as a tight-binding model of graphene the terms correspond to physically distinct gapping instabilities at the Dirac cones [37]. The mass term mS h of (1) corresponds to a charge density wave in which electrons preferentially sit on one of the two sub-lattices A or B [38], whereas a linear combination of m 3 S 3 (14) and m 5 S 5 (17) yields a bond density wave in which electrons are distributed on both A and B sublattices in a Kekulé texture [39].
Note that a mass term proportional to the bilinearΨγ 3 γ 5 Ψ, whether hermitian or antihermitian, is qualitatively different. In terms of two-component spinors the resulting action reads since all elements are proportional to the combinationD +m there is no parity operation realisable as a linear transformation on Ψ,Ψ which flips the sign ofD but leaves the mass term invariant. This corresponds to the non-time-reversal invariant "Haldane mass", realised in graphene-like systems by alternately circulating currents in adjacent half-unit cells [40].
In summary, there are three linearly independent, physically indistinguishable parityinvariant mass terms available for continuum four-component Dirac fermions in 2+1d. We will see that this furnishes a non-trivial test for lattice fermion formulations, such as the domain wall formulation presented in Sec. 4, in which the matrices γ 3 and γ 5 appear in inequivalent ways. Of course, the transformations (12,15) are merely special cases of the rotations (6,5) with α = π 4 , so the physical equivalence of the mass terms is a consequence of U(2) symmetry.

Ginsparg-Wilson Relations and the Overlap operator
The Nielsen-Ninomiya theorem [41,42] famously asserts the impossibility of a lattice fermion formulation with a physical (ie. undoubled) spectrum which simultaneously respects locality, unitarity (or reflection positivity in Euclidean metric) and the existence of a conserved axial charge, expressed in four dimensions via the anticommutator {D, γ 5 } = 0. Ginsparg and Wilson [30] proposed a minimal modification to these criteria for lattice chiral fermions to be viable, namely that chiral symmetry now be constrained by the GW relation where a is the lattice spacing. The non-vanishing right hand side of (19) is effectively an O(a) contact term for the anticommutator of the propagator D −1 with γ 5 , which should become unimportant in the long-wavelength limit aD → 0. The GW relation appropriate to odd-dimensional fermions in irreducible representations of the spinor algebra was investigated by Bietenholz and Nishimura [43], who found that the correct parity anomaly is recovered due to non-invariance of the fermion measure under a "generalised parity" transformation. For the reducible 2+1d theories of interest in this paper the GW relations generalise to: In addition we require the hermiticity relations The GW symmetries (19)(20)(21) are realised by the overlap operator [33] where the matrix A ≡ aD W ils − 1 1, where D W ils is the lattice Dirac operator for orthodox Wilson fermions. It has the properties It is manifest that γ 3 and γ 5 are treated identically in the overlap formalism. A distinct overlap operator, appropriate for irreducible fermions and reproducing the correct parity anomaly, was considered in [44].
Analogously to the situation in 4d [45] the GW relations (19-21) admit the following remnant symmetries, which recover the desired U(2N f ) symmetries (3)(4)(5)(6) in the continuum and long-wavelength limits a → 0, aD → 0: For theories such as the GN and Thirring models defined near a critical point it is not a priori clear that continuum and long-wavelength limits coincide.
To construct mass terms, we define projection operators of two kinds Relations (19)(20)(21) may be used to showγ 2 5 =γ 2 3 = 1, so that theP s are genuine projectors. Importantly, it also follows that Therefore it makes sense to define two different sets of projected fields: enabling a decomposition of the kinetic term in two different ways: The mass term in this approach is defined in terms of projected fields: Although this term only corresponds with the desired mass term in the continuum limit a → 0, the extra O(a) piece has the same form as the kinetic term, so the overall effect is a benign wavefunction renormalisation. We can also check the effects of the GW symmetries (24), along with the fermion number symmetry (3). With obvious notation, and to O(α): The variation (33) provides the interpolators for the Goldstone modes associated with a spontaneous breaking Ψ (1 − aD 2 )Ψ = 0; we see there is an O(a 2 ) correction to the expected continuum forms.
It is interesting to repeat this exercise for the two other parity-invariant mass terms in (14,17). The mass terms constructed from the projected fields are In this case the O(a) correction does not correspond to an existing term in the action, so the correction may be less benign. Another symptom is that if we label the full Dirac operator including mass terms as D = D + m h S 1 + m 3 S 3 + m 5 S 5 , then but that no such simple relations exist once m 3 , m 5 = 0. We find the O(α) variations of the mass term S 3 under (24) with similar results, mutatis mutandis for variations of S 5 . In all cases the expected expressions for Goldstone interpolators are recovered in the continuum limit. However for a > 0 there are subtle differences in the effect of the GW variations on the different mass terms; for instance, the two Goldstone interpolators "ψψ" and "ψγ 5 ψ" resulting from δ GW in (36) have realisations differing by O(a 2 ) away from the continuum limit. In this sense, full U(2) symmetry is only recovered by overlap fermions as a → 0.

Domain Wall Formulation
In this section we set out an alternative route to undoubled U(2N f )-symmetric fermions using the domain wall approach first introduced by Kaplan [31] and put into the present form by Furman and Shamir [32]. Here we follow the treatment set out in [17]. Define 4-spinor fields Ψ(x, s),Ψ(x, s) living on a 4d lattice where x denotes the usual 2+1d spacetime coordinates of a lattice site and s = 1, . . . , L s its coordinate along the extra dimension, here labelled 3. The kinetic term in the action is then x,y s,rΨ (x, s)D DW (x, s|y, r)Ψ(y, r).
with domain wall Dirac operator The first term is the orthodox 2 + 1d Wilson operator with gauge link variables U µ (x), and D DW 3 controls hopping in the 3 direction: Note there are Dirichlet boundary conditions imposed in direction 3, at s = 1 and s = L s . The inclusion of D DW 3 explicitly destroys the equivalence of γ 3 and γ 5 in the dynamics described by the action (37), so it will be important to test whether and how this is recovered in practice.
The key idea [31] is that the dynamics generated by (39) and (40), with suitably chosen M, results in fermion zeromodes localised on domain walls at s = 1, L s , which are also respectively ∓ eigenmodes of γ 3 . The 2 + 1d physics we wish to describe is formulated entirely using these localised modes (the Wilson terms in (39,40) render the would-be zeromodes due to unwanted doubler species non-normalisable in the limit L s → ∞ [31]). In particular we need to define 2 + 1d fermion mass terms corresponding to their continuum counterparts in (1), (14) and (17). To this end, define fermion fields ψ(x),ψ(x) living in 2+1d: where from now on P ± ≡ 1 2 (1 ± γ 3 ). We thus consider actions of the form (37) supplemented by three alternative mass terms: It is interesting to note that S h has the same form as the fermion mass term for domain wall formulations of 3+1d physics, and couples fields from opposite walls; S 3 also couples opposite walls, but S 5 couples fields living on the same wall.
In the next section we will examine the numerical consequences of the three terms (42)(43)(44) and in particular check whether they yield compatible, U(2)-symmetric results in the L s → ∞ limit.

Numerical Results
In order to explore the domain wall action (37) supplemented by one of the mass terms (42)(43)(44) we have performed quenched simulations of non-compact QED, so that U µ (x) = exp(iθ µ (x)) with the real link variable θ equilibrated using the Maxwell-like action The dimensionless coupling β ≡ (g 2 a) −1 , where the dimensionful fermion charge g is the natural scale with which to define physical observables. The continuum limit lies at β → ∞. In 2 + 1d QED is an asymptotically-free theory whose infra-red behaviour is still imperfectly understood. Since the non-compact lattice formulation is non-confining, numerical simulations are plagued by the slow fall-off of the photon propagator (∝ r −2 for the quenched theory), and the thermodynamic limit is extremely difficult to achieve at weak coupling. Quenched simulations with staggered lattice fermions presented in [35] suggest that chiral symmetry is spontaneously broken in the ground state in the continuum thermodynamic limit, namely (in lattice units a = 1) Lattice sizes up to 80 3 were used in support of this claim. Away from the continuum limit massless staggered fermions have a manifest U(1)⊗ U(1) global symmetry which spontaneously breaks to U(1) either when an explicit mass is introduced or a chiral condensate of the form (46) develops. Only in the continuum limit can this pattern enlarge to the expected U(2N f ) →U(N f )⊗U(N f ), with N f = 2 [20]. The main aim of this study is not to verify (46) for domain wall fermions, but rather to test the restoration of the correct symmetry breaking pattern, with N f = 1, as L s → ∞. To this end we have performed simulations on system sizes 24 3 × L s , with couplings β = 0.5, 1.0 and 2.0, corresponding respectively to strong, intermediate and weak coupling regimes. Throughout the domain wall fermion mass M in (39) was set to 0.9 (physical quantities are expected to be M-independent; reflection positivity requires 0 < M < 1), and unless otherwise stated the masses m h , m 3 , m 5 in (42-44) to 0.01. For convenience a hybrid Monte Carlo algorithm was used to generate the {θ} ensembles, with 100 trajectories of average length 1.0 between measurements; a much more efficient Fourier-space method was employed in the original study [35].

Bilinear Condensates
First we explore bilinear condensates generically defined by ψ Γ i ψ ≡ ∂ ln Z/∂m i , where m i ∈ {m h , m 3 , m 5 }. These form in response either to m i = 0 or in the thermodynamic limit to strong dynamics. For the latter case in the limit m i → 0 all three condensates should be equal if U(2) symmetry is manifest. It is not obvious that this will occur for the action (37), both because of the Wilson terms inherent in (39,40), and because γ 3 and γ 5 do not appear in (37) on an equal footing. By hypothesis, however, U(2) symmetry should be recovered as L s → ∞.
To begin, we present results obtained using a spatial point source on a configuation generated at β = 0.5 (in fact, the numbers result from averaging over 10 spatial sources); the systematics are easiest to expose at the strongest coupling. Note from (42)(43)(44) that each condensate gets contributions from two terms: for ψ ψ and i ψ γ 3 ψ the two terms arise from four-dimensional propagators running from s = 1 to L s and L s to 1 respectively; for i ψ γ 5 ψ each contribution is from a propagator starting and ending on the same domain wall. Within the working numerical precison each contribution is the complex conjugate of the other, so the sum is real. However, it turns out the imaginary component parametrises the approach to the U(2)-symmetric limit. Define i Ψ (1)γ 3 Ψ(L s ) = i 2 ψ γ 3 ψ Ls + i∆ h (L s ) (where the first term is real, and the spatial coordinate x is suppressed), and then write: The residuals ∆ h and ǫ i must vanish for a U(2)-invariant limit.  note that ∆ h is measured directly as the imaginary component of i ψ γ 3 ψ using just the + components of Ψ,Ψ, while to estimate the ǫ i the value of i ψ γ 3 ψ Ls→∞ is taken to be that measured at L s = 48. Several features are apparent: • The dominant correction by almost an order of magnitude is ∆ h , which contributes to the hermitian condensate ψ ψ but not, as a result of the twist, to the antihermitian i ψ γ 3,5 ψ . Indeed, at the weakest coupling β = 2.0 ∆ h is the only residual large enough to measure.
• All residuals decay approximately as exp(−cL s ), consistent with U(2) symmetry restoration as L s → ∞.
• ǫ 3 and ǫ 5 are practically identical -the disparity in Fig. 1 at L s = 40 probably results from uncertainty in the determination of i ψ γ 3 ψ L S →∞ . This is striking since the underlying structures in terms of four-dimensional propagators are quite distinct.
Of course, we should not draw universal conclusions from a single gauge configuration, particularly since the bilinear condensates display strong intermittent upward fluctuations across a gauge ensemble. Fig. 2 shows ∆ h (L s ) evaluated for 1200 (24 3 ) or 2400 (12 3 ) measurements of i ψ γ 3 ψ at the three couplings explored (for L s ≥ 16 O(500) measurements were made at the strongest coupling β = 0.5 on 24 3 ) using a stochastic noise vector η(x); to faithfully implement the projectors P ± in (43) the same η is chosen for all spin components in the trace. Reassuringly, for L s > ∼ 8 the exponential behaviour ∆ h ∝ e −cLs persists, with decay constant c decreasing as the coupling increases. Comparison of the β = 0.5 results also shows c is in general also volume-dependent. Although not investigated here, we also expect c to depend on the Lagrangian parameter M. There is no reason to doubt, however, that the residuals will eventually become an insignificant systematic as L s → ∞. Results for the condensate i ψ γ 3 ψ , which as a result of the twist does not include effects of the dominant residual ∆ h , are shown in Fig. 3. With the use of a stochastic estimator a conjugate gradient residual of 10 −12 per (four-dimensional) lattice site and spin component sufficed. As a function of L s the results plateau essentially once ∆ h /| ψ γ 3 ψ | is smaller than the statistical error confirming the trends of Fig. 2: this occurs for L s > ∼ 8 for β = 2.0 and L s > ∼ 16 for β = 1.0. At the strongest coupling β = 0.5 the available numerical resources have not permitted this regime to be probed on the reference 24 3 volume; however data taken on 12 3 eventually plateau for L s > ∼ 28. Based on the ∆ h data from Fig. 2 we might guess the β = 0.5, 24 3 plateau may set in for L s ≈ 40. Fig. 3 also shows results for the other bilinear condensates ψ ψ (h) and i ψ γ 5 ψ (5) for β = 1.0, showing that within statistical errors they become consistent with i ψ γ 3 ψ for L s > ∼ 20. The impact of the ∆ h residual on ψ ψ is clearly discernable, and suggests it is not the optimal choice for systematic studies of chiral symmetry breaking. However, with this level of precision there is no reason not to suppose U(2) symmetry is eventually restored. Fig. 4 shows the number of conjugate gradient iterations required for each fermion matrix inversion is roughly constant once the plateau is reached, implying that the numerical effort needed scales linearly with L s as the U(2)-symmetric limit is approached.
Finally, Fig. 5 plots i ψ γ 3 ψ as a function of m 3 at β = 1.0 on 24 3 × 20. The data can be plausibly extracted linearly to the chiral limit m 3 → 0 to yield a non-vanishing intercept, consistent with the hypothesis that chiral symmetry is indeed spontaneously broken (though at this stage no thermodynamic or continuum limit is taken). Also shown is the staggered fermion condensate 1 2 χχ (ie. normalised to N f = 1) on the same spatial volume estimated from Fig. 7 of [35]. The values are close enough to be plausibly consistent, but in any case provide reassurance first that the domain wall fermion data do not contain significant contributions from doubler species, and secondly that the explicit chiral symmetry-breaking violation by the Wilson terms in (39) is mitigated in the domain wall approach. Of course, since there are as-yet unquantified renormalisations of both elementary fields and composite operators, these considerations strictly require a continuum limit to be taken.

Meson Propagators
Next we turn attention to correlations between mesonic operatorsψΓψ at different spacetime points, with Γ ∈ {γ 5 , γ 3 , 1 1, iγ 3 γ 5 } all yielding spin-0 states, two with positive parity and two negative. We again focus on the recovery of U(2) symmetry as L s → ∞. A similar analysis for staggered fermions in the continuum limit of non-compact QED 3 was presented in [46]. We begin by listing two important identities of the 2+1+1d fermion propagator S(m i ; x, s; y, s ′ ) = Ψ(x, s)Ψ(y, s ′ ) , which follow directly from (37) and have been checked numerically on small lattices: Here x denotes the coordinate in 2+1d and s that in the 3rd direction,s ≡ L s −s+1, and the dagger acts on Dirac indices. In what follows a mass parameter m i is omitted from the argument list if it takes the value zero. The identities (50,51) are of course crucial for the efficient calculation of timeslice correlators with a minimal number of matrix inversions, although there will be instances where a further inversion with the sign of the mass term m 3,5 reversed is needed. The contrast with relations (35) for overlap fermions, which away from the continuum limit only hold for m 3 = m 5 = 0 should be noted. As will be demonstrated, in the limit L s → ∞ two further approximate relations become valid: and equivalent ones for 3 ↔ 5. Now consider the pion correlator C 5 (x) = ψ (0)γ 5 ψψ(x)γ 5 ψ(x) . Under the parity definition (10) this state has J P = 0 − . Using the definition (41) and ignoring the overall minus sign from the Grassmann nature of Ψ,Ψ, we deduce the following relation in terms of 2+1+1d propagators C 5 (x) = tr P + S(m i ; x, L s ; 0, L s )P − γ 5 P + S(m i ; 0, L s ; x, L s )P − γ 5 + P − S(m i ; x, 1; 0, 1)P + γ 5 P − S(m i ; 0, 1; x, 1)P + γ 5 + P − S(m i ; x, 1; 0, L s )P − γ 5 P + S(m i ; 0, L s ; x, 1)P + γ 5 It's natural in this formalism to have projectors P ± sandwiching the fermion propagators -they are the bridge, both formally and in the code, between the 4d and 3d worlds. Now specialise to the case of mass term S h , and use the identity (50), along with translational invariance in 3d, to arrive at C h 5 (x) = tr S(m h ; 0, L s ; x, L s )P − S † (m h ; 0, L s ; x, L s )P + + S(m h ; 0, 1; x, 1)P + S † (m h ; 0, 1; x, 1)P − + S(m h ; 0, 1; x, L s )P − S † (m h ; 0, 1; x, L s )P − + S(m h ; 0, L s ; x, 1)P + S † (m h ; 0, L s ; x, 1)P + (54) This has the familiar look of a pion correlator (ie. positive definite terms of the form SS † ) apart from the projectors, and is clearly calculable with two sources, one located at s = 1 and the other at s = L s . In the abbreviated form (55) note that −+ and +− correlators involve fermion propagators linking a domain wall to itself, whereas for −− and ++ the propagators run between the walls. Next consider the 0 + correlator C h 3 (x) = ψ (0)γ 3 ψψ(x)γ 3 ψ(x) . This produces an expression with a different 4d spacetime structure, but which with the help of (51) can be rearranged to coincide exactly with (54): These two states are thus exactly degenerate if the hermitian mass term S h is chosen. The same methods can be used for the correlators which are respectively 0 + and 0 − to show Therefore these two states are also degenerate, but distinct fromψγ 3 ψ andψγ 5 ψ in (54). Assuming the primitive correlators C h±± are all positive definite, then |C h If the symmetry breaking in the direction of S 3 is spontaneous this is consistent with U(2)→U(1)⊗U(1) yielding two light pseudo-Goldstone modes interpolated byψγ 5 ψ andψγ 3 ψ. Had we instead chosen (11) as the definition of parity, then the assignments of the two Goldstone states would be reversed so thatψγ 5 ψ is now 0 + andψγ 3 ψ 0 − . However, the overall picture of a Goldstone 0 ± doublet and a non-Goldstone 0 ± doublet remains unaltered.
The picture is more interesting when the mass term S 3 is considered. The correlator C 3 5 = C 3−+ + C 3+− + C 3−− + C 3++ in analogy to (54,55), but for C 3 3 we find The expression (60) is similar in form to (56) but requires twice the number of matrix inversions to evaluate. However, use of the large-L s approximations (52) gives Similarly we find This suggests that the Goldstone modes are now interpolated byψψ (0 + ) andψγ 5 ψ (0 − ), consistent with the variation of S 3 under (3-6), but that the U(2)→U(1)⊗U(1) pattern will now only be recovered in the limit L s → ∞.
Finally, for mass term S 5 we have the apparent contradiction with the expected identification ofψψ (0 + ) andψγ 3 ψ (0 − ) (using either the symmetries (3)(4)(5)(6) or the change of variables (15)) as Goldstone interpolators will be further discussed below. We now present numerical results for the primitive correlators C ±± (x). It is clear the recovery of U(2) symmetry hinges on the validity of the approximation (52). As in Sec. 5.1, the first test uses the numerically most demanding system of β = 0.5 on 24 3 × L s , with fermion mass m i = 0.01. For a single equilibrated configuration we  Fig. 6 for mass terms m 3 S 3 , m h S h , m 5 S 5 , along with the same quantity evaluated for |C 3 | (labeled "3-" in the plot). It is apparent that C 3±− has reached its large-L s limit by L s ≈ 40, but that C h±− converges rather more slowly. Most importantly, C 3+− is indistinguishable from C 3+− , whereasC 3−− actually changes sign around L s ≈ 25 and by L s ≈ 40 is practically equal to −C 3−− , consistent with (52). The correlators C 5±− do not, however, fit the same pattern; C 5+− lies systematically above C 3+− , while C 5−− → 0 as L s → ∞.  This picture is confirmed by a study of 500 configurations at β = 1.0, m i = 0.01 on 24 3 . Only the components C ±− , which are manifestly real, are calculated; the counterparts C ∓+ are identical, and the correlatorsC ±∓ requiring an inversion with a negative fermion mass recover these properties as L s → ∞. To focus on systematic effects the same ensemble was used for each L s -the resulting statistical errorbars show the expected evidence of autocorrelation between timeslices. Fig. 7 compares C 3+− with C h+− ; C 3+− converges to its large-L s limit for L s > ∼ 12, whereas C h+− converges somewhat more slowly, requiring L s > ∼ 16. Crucially though, by L s = 24 the two are indistinguishable. Fig. 8 confirms that after a sign change occuring for L s ≈ 8 then lim Ls→∞C 3−− = −C 3−− , and that both also coincide in magnitude with C h−− for L s > ∼ 24. We have thus demonstrated that the mass terms S h via relations (56,58), and S 3 via (60-63), must yield equivalent meson spectra as L s → ∞. Finally Fig. 9 shows the primitive correlators C 5±− (t) for various L s . This time the picture is different: lim Ls→∞ C 5+− = C 3+− + C 3−− , while lim Ls→∞ C 5−− = 0. The relations (64-67) then predict that all four meson states become degenerate in this limit, and in this sense U(2) symmetry is realised. Explicit breaking to U(1)⊗U(1) by a mass term m 5 S 5 is not accomplished, however. We can trace this back to the expressions for the would-be Goldstones interpolated byψγ 3 ψ (65) andψψ (66) which are exact even for finite L s , and can only coincide in magnitude in the limit that either C 5−− or C 5+− vanishes. There is thus no opportunity for degeneracy breaking induced by S 5 to develop; the system responds to this constraint by C 5+− becoming larger as L s → 0 so that all four mesons become degenerate with the two would-be Goldstones observed with mass terms m h S h , m 3 S 3 . Whether this exceptional behaviour is restricted to a singular point on the U(2) manifold or whether it is a more general phenomenon must be the subject of further study.

Discussion
The main focus of this paper is the numerical investigation of domain wall fermions for problems involving fermions in reducible spinor representations in 2+1d. The results of Sec. 5 suggest it is highly plausible that the correct global U(2) symmetries are recovered in the limit L s → ∞, independently of the continuum limit, just as for QCD [32]. How quickly the limit is approached remains, of course, a dynamical issue. However, the enhanced global symmetries of the problem admit new, helpful features not present in QCD, namely the ability to introduce a U(2)-breaking mass term in any of three independent directions all yielding equivalent results as L s → ∞, aiding the robust determination of the L s required in practice. It is also significant that the dominant artifact ∆ h (L s ), related to the "residual mass" in QCD simulations, cancels from the physical observables explored here when the twisted form im 3,5ψ γ 3,5 ψ is used.
The only exception to the general good news is the failure of the meson spectrum calculated with im 5ψ γ 5 ψ to manifest the expected U(2)→U(1)⊗U(1) pattern, essentially because the formalism results in distinct expressions for the Goldstone correlators in this (possibly singular) case. The system responds by forcing their difference to vanish as L s → ∞, so that no symmetry breakdown occurs. The general conclusion, however, must be that the method looks extremely promising, and optimal if symmetry breaking is implemented with a mass term im 3ψ γ 3 ψ. The next step is to explore a fully dynamical implementation of domain wall fermions, and then test U(2) restoration by computing the separate components of the axial Ward identities. In this regard it is worth noting that the identity (50) implies that the fermion determinant is real for mass terms S h and S 3 , and hence there is no Sign Problem obstruction to running hybrid Monte Carlo with even N f . For mass term S 5 the identity (51) does not suffice to prove detD DW real, but is completely analogous to the corresponding identity for domain wall fermions in 4d.
An interesting distinction with QCD, or with more general 4d theories, is the relation between GW/overlap and domain wall approaches. As stressed, the two approaches differ in that GW formulation presented here maintains strict equivalence between γ 3 and γ 5 whereas the domain wall treats them very differently. We have seen in Sec. 3 that for GW fermions the full U(2) symmetry is only recovered in the continuum limit: this is seen either via the effect of the remnant symmetry rotations δ GW on bilinears (36), or from the failure of the fermion propagator identity (35) to generalise to m 3,5 = 0. Perhaps the problem arises from the GW implementation of the twisted mass terms (34), which as noted introduces an O(a) correction which cannot be compensated by a simple field rescaling. On the face of it, then, GW/overlap fermions can only recover U(2) for a → 0, whereas domain wall fermions apparently recover U(2) for L s → ∞ irrespective of whether a continuum limit is taken. One might speculate that there is a different set of GW relations distinct to (19,20,21), in which γ 3 and γ 5 do not appear on an equivalent footing, but which is consistent with the large-L s limit of domain wall fermions even away from the continuum limit. It would also be very valuable to understand things from a more fundamental perspective by attempting to derive the overlap from the large-L s limit of domain wall fermions, as done for 4d in [47]. It will also be interesting to extend consideration to the case of non-vanishing chemical potential µ, where there has been some discussion over the correct implementation of the GW symmetries [48,49]; happily, there are several interesting 2+1d systems where simulations at µ = 0 with orthodox Monte Carlo methods are feasible [50,51,52].
The biggest goal, however, once dynamical fermions have been implemented, is to study the strongly coupled fixed points of the GN, Thirring and graphene [53,54] systems to see if the critical exponents match those previously obtained with staggered fermions, and to check whether the expected dependence on N f is seen (an exploratory study of the GN model with domain wall fermions is presented in [55]). This would not only provide an important and hitherto unexplored test of the applicability of the domain wall approach to lattice fermions, but more generally would represent an important step in the definition of interacting quantum field theories beyond weak coupling.