Numerical study of the $2+1d$ Thirring model with U($2N$)-invariant fermions
Simon Hands










Department of Physics, College of Science,

Swansea University, Singleton Park, Swansea SA2 8PP, United Kingdom


In 2+1 dimensions the global U(2N) symmetry associated with massless Dirac fermions is
broken to U(N)⊗U(N) by a parity-invariant mass. I will show how to adapt the
domain wall formulation to recover the U(2N)-invariant limit in interacting
fermion models as the domain wall separation is increased. In particular, I will
focus on the issue of potential dynamical mass generation in the Thirring model,
postulated to take place for N less than some critical N_c. I will present
results of simulations of the model using both HMC (N=2) and RHMC (N=1)
algorithms, and show that the outcome is very different from previous numerical
studies of the model made with staggered fermions, where the corresponding
pattern of symmetry breaking is distinct.



Numerical study of the 2+1d Thirring model
 with U(2N)-invariant fermions

    
Simon Hands1Acknowledges financial support
from  STFC and the Leverhulme Trust

    
====================================================================================




§ INTRODUCTION


Relativistic fermions moving in 2+1 spacetime dimensions form
the basis of some facscinating quantum field theories, whose principal 
applications mainly fall in the domain of the condensed matter physics of layered
systems. For example, fermion degrees of freedom described by a
Dirac equation have been invoked and studied in models of the underdoped regime
of d-wave superconductors <cit.>, 
where they arise as excitations at the nodes of the
gap function Δ(k⃗); in models of spin-liquid behaviour in Heisenberg
antiferromagnets <cit.>, 
where the infrared behaviour of compact QED_3 remains an
open issue; as surface states of three-dimensional topological insulators; and of
course as low-energy electronic excitations in graphene <cit.>.

What all these models share in common is their use of four-component reducible spinor representations for the fermi fields. For non-interacting 
fermions the action in Euclidian metric is 

    S=∫ d^3x ψ̅(γ_μ∂_μ)ψ+mψ̅ψ;  μ=0,1,2;

a key point is the mass term proportional to m is hermitian and invariant
under the parity inversion x_μ↦-x_μ. For m=0
S has a global U(2) invariance generated by 

    ψ↦ e^iαψ, ψ̅↦ψ̅e^-iα;      ψ↦ e^αγ_3γ_5ψ, ψ̅↦ψ̅e^-αγ_3γ_5;


    ψ↦ e^iαγ_3ψ, ψ̅↦ψ̅e^iαγ_3;   ψ↦ e^iαγ_5ψ, ψ̅↦ψ̅e^iαγ_5.

For m≠0 the γ_3 and γ_5 rotations (<ref>) are no
longer symmetries and the general pattern of breaking is thus
U(2N)→U(N)⊗U(N), where we generalise to N degenerate flavors.

Because there is no chiral anomaly in 2+1d, it is possible to perform a change
of variables in the path integral to identify two further “twisted” mass terms which,
though antihermitian, are physically equivalent to the mψ̅ψ  of (<ref>):

    im_3ψ̅γ_3ψ;   im_5ψ̅γ_5ψ.

The “Haldane” mass term m_35ψ̅γ_3γ_5ψ is not
equivalent because it changes sign under parity, and will not be considered
further.





§ THE THIRRING MODEL IN 2+1D


A model of particular interest is the Thirring model, which has a contact interaction
between conserved fermion currents. Its Lagrangian density reads

    L=ψ̅_i(∂/ +m)ψ_i+g^22N(ψ̅_iγ_μψ_i)^2;  i=1,…,N.

Equivalently, its dynamics is captured by a bosonised action via the introduction
of an auxiliary vector field A_μ:

    L_A=ψ̅_i(∂/ +ig√
    NA_μγ_μ+m)ψ_i+12A_μ^2.

The Thirring model is arguably the simplest interacting theory of fermions
requiring a computational approach. The coupling g^2 has dimension 2-d,
where d is the dimension of spacetime, so a naive expansion in powers of g^2
is non-renormalisable for d>2. However, things look different after
resummation. First introduce an additional Stückelberg scalar φ so
the bosonic term becomes 12(A_μ-∂_μφ)^2, to
identify a hidden local symmetry <cit.>

    ψ↦ e^iαψ;  A_μ↦
    A_μ+∂_μα;  φ↦φ+α.

This point of view strongly suggests the identification of A_μ as an abelian
gauge field; the original Thirring model (<ref>) corresponds to a
unitary gauge φ=0. In Feynman gauge, to leading order in 1/N the
resummed vector propagator is now of the form ⟨
A_μ(k)A_ν(-k)⟩∝δ_μν/k^d-2. An expansion in powers
of 1/N may now be developed by analogy with QED_3, and is exactly
renormalisable for 2<d<4 <cit.>.

The outstanding theoretical issue is whether, for g^2 sufficiently large and
N sufficiently small, there is a symmetry-breaking transition leading to
formation of a bilinear condensate ⟨ψ̅ψ⟩≠0
accompanied by dynamical fermion mass generation. 
It is of particular
interest to identify the critical N_c above which symmetry breaking does not
occur even in the strong coupling limit. 
The transitions at
g_c^2(N<N_c) then potentially define a series of distinct quantum critical
points (QCPs). One early prediction, using strong-coupling Schwinger-Dyson
equations in the ladder approximation, found N_c≃4.32 <cit.>.



§ THE THIRRING MODEL WITH STAGGERED FERMIONS


The Thirring model has been studied by numerical simulations using staggered
fermions, with action <cit.>

    S_stag=12∑_xμ i[χ̅_x^iη_μ x(1+iA_μ
    x)χ_x+μ̂^i -
    ]+m∑_xiχ̅_x^iχ_x^i+N4g^2A_μ x^2.

Eq. (<ref>) is not unique, but has the feature that the linear coupling
of the auxiliary precludes higher-point interactions between
fermions once it is integrated over. The action (<ref>) has a global
U(N)⊗U(N) symmetry broken to U(N) by either explicit
or spontaneous mass generation. In a weakly coupled long-wavelength limit a
U(2N_f) symmetry is recovered with N_f=2N <cit.>; however the putative QCP is not
weakly coupled, so this conventional wisdom must be questioned. 






The action (<ref>) has been studied for
N_f∈[2,6] <cit.> and QCPs have
been identified and characterised. Other studies with N_f=2 have found
compatible results <cit.>. 
A more recent study which plausibly identifies the
strong-coupling limit reported N_fc=6.6(1), δ(N_fc)≈7 <cit.>. 
The results are summarised in
Fig. <ref>. The “chiral” symmetry is indeed broken for small N_f and large
g^2, and the exponent δ defined by the critical scaling
⟨χ̅χ⟩|_g_c^2∝ m^1δ is very
sensitive to N_f. Other exponents can be estimated using hyperscaling.
These results are in qualitative but not
quantitative agreement with the Schwinger-Dyson approach, which predicts
δ(N_fc)=1. 
A non-covariant form of (<ref>) has been used to model the
semimetal-insulator transition in graphene, finding
N_fc≈5 <cit.>
suggesting that a Mott insulating phase is
possible for the N_f=2 appropriate for monolayer
graphene <cit.>.

The staggered Thirring model thus exhibits a non-trivial phase diagram, with a
sequence of QCPs with an N-dependence  
quite distinct from those of the
theoretically much better-understood Gross-Neveu (GN) model. However, recent
simulations of N_f=2 with a fermion bag algorithm which permits study directly
in the massless limit, have found compatible exponents for the 
QCP <cit.>:

    ν   =   0.85(1);     η=0.65(1);     η_ψ=0.37(1);     
    ν   =   0.849(8);   η=0.633(8);   η_ψ=0.373(3);

The large-N GN values are ν=η=1. These results are troubling: from
the perspective of the large-N expansion using bosonised actions the models
should be distinct, whereas (<ref>) suggest rather they lie in the
same RG basin of attraction. Indeed, when written purely in
terms of four-point interactions between staggered fields spread over elementary
cubes, the only difference between the models is an extra body-diagonal coupling
in the GN case <cit.>.  

In this study we examine the possibility that staggered
fermions do not reproduce the expected physics because of a failure to
capture the correct continuum global symmetries near a QCP. A similar insight has
been offered by the Jena group <cit.>. 




§ DOMAIN WALL FERMIONS IN 2+1D


The physical idea of domain wall fermions (DWF) is that fermions Ψ,Ψ̅ are allowed to
propagate along an extra fictitious dimension of extent L_s with open boundary
conditions. In 2+1+1d this propagation is governed by an operator
∼∂_3γ_3. As L_s→∞ zero modes of D_DWF localised on the domain walls at
either end become ± eigenmodes of γ_3, and physical fields in the
target 2+1d space identified via

    ψ(x)=P_-Ψ(x,1)+P_+Ψ(x,L_s);   ψ̅(x)=Ψ̅(x,L_s)P_-+Ψ̅(x,1)P_+,

with P_±=12(1±γ_3). The walls are then coupled with a term
proportional to the explicit massgap m. 

However, the emergence of the U(2N) symmetry outlined in Sec. <ref> is
not manifest, because while the wall modes are eigenmodes of γ_3, the
continuum symmetry (<ref>,<ref>) demands equivalence under
rotations generated by both γ_3 and γ_5. Another way of seeing
this is that the twisted mass terms (<ref>) should yield identical
physics, eg. the strength of the corresponding bilinear condensate, as
L_s→∞. This requirement is apparently non-trivial, since while m and
m_3 couple Ψ and Ψ̅ fields on opposite walls, m_5 couples
fields on the same wall.

The recovery of U(2N) symmetry as L_s→∞ was demonstrated numerically
in a study of quenched non-compact QED_3 on 24^3× L_s systems, for
a range of couplings β <cit.>. First define the principal residual Δ
via the imaginary part of the twisted condensate:

    [i⟨Ψ̅(1)γ_3Ψ(L_s)⟩]=
    -[i⟨Ψ̅(L_s)γ_3Ψ(1)⟩]≡Δ(L_s).

The difference between the various condensates and their value
in the large-L_s limit is then specified in terms of secondary residuals
ε_i(L_s) via

    ⟨ψ̅ψ⟩_L_s   =   i⟨ψ̅γ_3ψ⟩_L_s→∞+2Δ(L_s)
    +2ε_h(L_s);
    
    i⟨ψ̅γ_3ψ⟩_L_s   =   i⟨ψ̅γ_3ψ⟩_L_s→∞
    +2ε_3(L_s);
    
    i⟨ψ̅γ_5ψ⟩_L_s   =   i⟨ψ̅γ_3ψ⟩_L_s→∞
    +2ε_5(L_s).


Empirically, as shown in Fig. <ref>, the residuals Δ and ε_i
decay exponentially with L_s, with a clear hierarchy
Δ≫ε_h≫ε_3≡ε_5. This is strong
evidence for the ultimate recovery of U(2N) (the convergence rate is strongly
dependent on both β and system volume), and moreover suggests the optimal
simulation strategy is to focus on the twisted condensate
i⟨ψ̅γ_3ψ⟩ for which  finite-L_s corrections are
minimal. The equivalence of γ_3 and γ_5 condensates at finite
L_s and their
superior convergence to the large-L_s limit was shown analytically in
<cit.>, where the convergence to fermions obeying 2+1d Ginsparg-Wilson
relations was also demonstrated. Exponential improvement of convergence was 
shown in the large-N limit of the GN model in <cit.>.  The benefits of a 
twisted mass term for improved recovery of U(2N) symmetry were also observed in
a study of non-compact QED_3 using Wilson fermions in <cit.>.





§ THE THIRRING MODEL WITH DOMAIN WALL FERMIONS


Even after settling on DWF, we still encounter some remaining
formulational issues. Just as for the staggered model, we have chosen a linear
interaction between the fermion current and a 2+1d vector auxiliary field
A_μ defined on the lattice links. The simplest approach to formulating the Thirring model is
to restrict the interaction to the physical fields (<ref>) defined on the domain walls at
s=1,L_s. This follows the treatment of the GN model with DWF developed in
<cit.>. It has the technical advantage that the Pauli-Villars fields
required to cancel bulk mode contributions to the fermion
determinant do not couple to A_μ, and hence can be safely excluded from the
simulation, which brings significant cost savings. In what follows this approach
will be referred to as the Surface model.

However, following the discussion below eq. (<ref>) suggesting the strong
similarity of A_μ with an abelian gauge field, we also consider a Bulk
formulation in which Ψ, Ψ̅ interact with a “static”
field (ie. ∂_3 A_μ=0) throughout the bulk:

    S=Ψ̅ DΨ=Ψ̅D_WΨ+Ψ̅D_3Ψ+m_iS_i,

with 

    D_W=γ_μ D_μ-(D̂^2+M);   D_3=γ_3∂_3-∂̂_3^2;

and m_iS_i is the explicit mass term defined only on the walls. Here 

    D_μ xy   =   12[(1+iA_μ x)δ_x+μ̂,y-(1-iA_μ
    x-μ̂)δ_x-μ̂,y],
    D̂^2_xy   =   12∑_μ[(1+iA_μ x)δ_x+μ̂,y+(1-iA_μ
    x-μ̂)δ_x-μ̂,y-2δ_xy],
  
and M is the domain wall height, here set equal to 1. ∂_3, ∂̂_3^2 are defined
similarly using finite differences in the 3 direction, respecting the open
boundary conditions, and with no coupling to
the auxiliary. These definitions imply the following properties:

    [∂_3,D_μ]=[∂_3,D̂^2]=0      
    [∂_3,∂̂_3^2]≠0

The action (<ref>) may be simulated using the HMC algorithm for N
even; however the failure of the third commutator in (<ref>) to vanish
everywhere is an
obstruction  to proving D is positive definite, so N=1 is
simulated using the RHMC algorithm with functional measure (
D^† D)^12.



§ NUMERICAL RESULTS


We have studied both surface and bulk formulations of the Thirring model in the
coupling range ag^-2∈[0.1,1.0] with first N=2 <cit.> and now N=1. The RHMC
algorithm used 25 partial fractions to estimate fractional powers of D^† D. Most results are obtained on 12^3×16 (N=2) or
12^3×8 (N=1). An exploration of volume and finite-L_s effects for
N=2 was presented in <cit.>. Summary results for the bilinear condensate
i⟨ψ̅γ_3ψ⟩ with m_3a=0.01 are shown in
Fig. <ref>. The bulk model shows
signifcantly enhanced pairing for g^-20.5 compared to the
surface model, and as might be anticipated pairing is greater for N=1
than N=2. This trend continues until a maximum at ag^-2≈0.2. In
previous work with staggered fermions this has been identified with the effective location of the
continuum strong coupling limit <cit.>, 
since for stronger lattice couplings there is a breakdown
of reflection positivity <cit.>. We are thus confident that the range of couplings
explored includes the strong coupling limit. 

Fig. <ref> shows the auxiliary action over the same coupling range, compared with the
free field value d2. The difference between surface and bulk models
is apparently very striking, but it should be borne in mind that this is really
a comparison of UV properties of two different regularisations of 
ostensibly the same continuum field theory.
Once again, there is tentative evidence for a change of behaviour of the N=1 bulk
model at ag^-2≈0.5.


Fig. <ref> shows the ratio
i⟨ψ̅γ_3ψ⟩/m_3χ_π 
obtained for N=2, with the transverse susceptibility 
defined 

    χ_π=N∑_x⟨ψ̅γ_5ψ(0)ψ̅γ_5ψ(x)⟩.

For a theory where the ψ,ψ̅ fields respect U(2N) symmetry, the
2+1d generalisation of the axial Ward identity predicts the ratio
to be unity. Fig. 5 shows that this requirement is far from being met, and that
further work is needed to understand and calibrate the identification of the
physical fields via relations such as (<ref>). Again, the disparity
between bulk and surface models is striking, with neither being obviously
preferred.


Finally we consider whether U(2N) is spontaneously broken at strong coupling.
In <cit.>
the bilinear condensate for N=2 was examined as a function of bare mass m
across a range of couplings, and in every case a linear scaling
⟨ψ̅ψ⟩∝ m was found, indicative of U(2N) symmetry
being manifest as m→0. We now extend this study to N=1.
Figs. <ref>,<ref> show 
⟨ψ̅ψ(m)⟩ for ma=0.01,…,0.05; with a trivial
rescaling implemented by choosing the abscissa as
m/⟨ψ̅ψ(am=0.01)⟩, data from the entire range of couplings
studied collapses onto a near linear curve, for both bulk (left) and surface
(right) models. There is no evidence for any singular behaviour associated with
a symmetry-breaking phase transition, and it seems safe to conclude
lim_m→0⟨ψ̅ψ⟩=0. Assuming this picture persists in the
large volume and L_s→∞ limits, this provides strong evidence that the
critical flavor number in the U(2N)-symmetric Thirring model is constrained by

    N_c<1.




§ SUMMARY AND OUTLOOK


It has been shown that it is feasible to use DWF to study U(2N)-symmetric
fermions in 2+1d, and that use of a twisted mass term ∼
im_3ψ̅γ_3ψ optimises the recovery of the symmetry as
L_s→∞. A study of the Thirring model at strong fermion self-couplings
then shows that DWF capture a very different physics to that described by
staggered fermions, which are governed by a different global symmetry away from
the weak-coupling long-wavelength limit. While it is still not possible to settle on the preferred
formulation of the strongly-coupled Thirring model, both the bulk and surface
versions presented here are in agreement that the critical flavor number
N_c<1. Fortunately this is compatible with the results obtained using a
distinct U(2N)-symmetric approach involving the SLAC
derivative <cit.>, also
presented at this conference <cit.>. Recent studies of U(2N)-symmetric QED_3,
an asymptotically-free theory, have also concluded
N_c<1 <cit.>. 
It is noteworthy that a disparity between DWF
and staggered fermions in a very different physical context, namely near a conformal
fixed point in a 3+1d non-abelian gauge theory, has also been reported <cit.>.

In future work it will obviously necessary to study the effects of
L_s→∞, V→∞, and of varying the domain wall height M. It will
also be valuable to study the locality properies of the corresponding 2+1d
overlap operator, furnishing a non-trivial test of the DWF approach in a new
physical context. Another question to ponder is what does chiral symmetry
breaking actually look like for 2+1+1d DWF? We are currently investigating
this issue by quenching the Thirring model. Finally, functional renormalisation
group studies indicate that in the hunt for a QCP it may be interesting to
include a U(2N)-invariant Haldane interaction 
-g̃^2(ψ̅γ_3γ_5ψ)^2 <cit.>; the control offered by DWF will make
this a straightforward exercise.






