Critical Behaviour in the Single Flavor Thirring Model in 2+1d
SSeptember 2020
Critical Behaviour in the Single FlavorThirring Model in 2+1 d Simon Hands a , Michele Mesiti b , and Jude Worthy aa Department of Physics, College of Science, Swansea University,Singleton Park, Swansea SA2 8PP, United Kingdom. b Swansea Academy for Advanced Computing, Swansea University,Bay Campus, Swansea SA1 8EN, United Kingdom.
Abstract:
Results of a lattice field theory simulation of the single-flavor Thirring modelin 2+1 spacetime dimensions are presented. The lattice model is formulated usingdomain wall fermions as a means to recover the correct U(2) symmetries of the continuummodel in the limit where wall separation L s → ∞ . Simulations on 12 , × L s , varyingself-interaction strength g and bare mass m , are performed with L s = 8 , . . . ,
48, and theresults for the bilinear condensate (cid:104) ¯ ψψ (cid:105) fitted to a model equation of state assuming aU(2) → U(1) ⊗ U(1) symmetry-breaking phase transition at a critical g c . First estimatesfor g − c a and critical exponents are presented, showing small but significant departuresfrom mean-field values. The results confirm that a symmetry-breaking transition doesexist and therefore the critical number of flavors for the Thirring model N c >
1. Resultsfor both condensate and associated susceptibility are also obtained in the broken phaseon 16 ×
48, suggesting that here the L s → ∞ extrapolation is not yet under control.We also present results obtained with the associated 2+1 d truncated overlap operator D OL demonstrating exponential localisation, a necessary condition for the recovery ofU(2) global symmetry, but that recovery of the Ginsparg-Wilson condition as L s → ∞ is extremely slow in the broken phase.Keywords: Lattice Gauge Field Theories, Field Theories in Lower Dimensions, GlobalSymmetries 1 a r X i v : . [ h e p - l a t ] N ov Introduction
The Thirring model is a quantum field theory of relativistic fermions interacting viaa contact between conserved currents. In this paper we will examine the model in2+1 d spacetime dimensions, and therefore further specify the use of reducible (ie. four-component) spinor fields, permitting the formulation of a parity-invariant mass term.The Lagrangian density reads L = ¯ ψ i ( ∂/ + m ) ψ i + g N ( ¯ ψ i γ µ ψ i ) (1)where µ = 0 , , i indexing N fermion species is implied. For a reduciblerepresentation of the Dirac algebra, there are two independent matrices γ and γ whichanticommute with the kinetic term in (1), which results in a U(2 N ) global symmetrygenerated by the set { , γ , γ , iγ γ } . For m (cid:54) = 0 this breaks as U(2 N ) → U( N ) ⊗ U( N ).The symmetry can also break spontaneously through generation of a bilinear conden-sate (cid:104) ¯ ψψ (cid:105) (cid:54) = 0; while there are clear analogies with “chiral” symmetry breaking, thisnomenclature should be eschewed in odd spacetime dimension.Spontaneous U(2 N ) symmetry breaking is believed to be theoretically possible forsufficiently large self-coupling g and sufficiently small N . Previous investigations ofthis question have used truncated Schwinger-Dyson equations [1, 2, 3], and FunctionalRenormalisation Group [4, 5, 6]. The prototype scenario has a symmetry breakingtransition at g c ( N ), where g c is a increasing function of N . A UV-stable renormalisationgroup fixed point can be defined as g → g c , where we identify a Quantum Critical Point (QCP) such that there exists an interacting continuum field theory solely specified bythe field content, dimensionality and pattern of symmetry breaking. The critical flavornumber N c required for symmetry breaking defined by g c ( N ) | N = N c → ∞ , which inprinciple need not be integer, is an important property of the model. Since symmetrybreaking 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-perturbativequantum field theory.The model also provides a natural arena for the application of lattice field theorymethods, and may well be the simplest fermion field theory requiring a computationalsolution. It turns out the choice of fermion discretisation has undue influence. Manyearly studies [7, 8, 9, 10, 11, 12] used staggered fermions, which naturally support asymmetry breaking U( N ) ⊗ U( N ) → U( N ) [13]. The results of these studies are broadlyconsistent; the most wide-ranging of them, exploiting plausible assumptions about the g → ∞ limit on a lattice, found N c = 6 . N < N c appear rather sensitive to N , suggesting the existence of a rich family ofdistinct QCPs.There are reasons to question whether this prediction for N c is correct; issues arisingfrom a lattice perspective are reviewed in [14]. Here we present a non-rigorous plausibilityargument for caution based on symmetry. The Thirring model in 2+1 d may be studiedanalytically in the limit of large N , and the mass M V of a vector f ¯ f bound state2nterpolated by the current density ¯ ψγ µ ψ is predicted in this limit to be [15] M V m = (cid:114) πmg ; lim g →∞ M V m = 0 . (2)Hence in the strong-coupling limit the Thirring model is a theory of conserved currentsinteracting via exchange of a massless spin-1 boson. Asymptotically the boson propa-gator switches from the canonical D µν ( k ) ∝ k − to ∝ k − ; this UV behaviour is exactlythat predicted for the photon of QED in the IR limit. Hence the Thirring QCP, if it ex-ists, could well be identical with the IR fixed point of QED , and the critical N c for boththeories should then coincide. The parallels between the Thirring model and abeliangauge theory are more apparent still once an auxiliary vector field A µ is introduced, asreviewed 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 , based on the conjecture f IR ≤ f UV , where f IR = − π lim T → F T ; f UV = − π lim T →∞ F T , (3)and F is the thermodynamic free energy density. For asymptotically-free theories suchas QED f UV is related to the count of non-interacting constituents: f UV = 34 × N + 1 . (4)Here is the appropriate factor for Fermi-Dirac statistics in 2+1 d , 4 is the number ofspinor components per flavor and 1 counts the single physical polarisation state of thephoton, which remains unconfined in QED . For a phase with spontaneously-brokensymmetry the number of weakly-interacting degrees of freedom is the Goldstone countplus the photon. Hence f IR = (cid:40) N + 1 U(2 N ) → U( N ) ⊗ U( N ); N + 1 U( N ) ⊗ U( N ) → U( N ) . (5)For QED , and by extension the continuum Thirring model, the conjecture thereforepredicts N c ≤ , whereas for the symmetry breaking dictated by the staggered formula-tion the equivalent bound is the much less stringent N c ≤
12. This disparity is a strongmotivation for exploring lattice fermions with the correct global symmetry.Recently this programme has been developed in two distinct directions. Latticefermions employing the SLAC derivative, which is non-local but manifestly U(2 N )-symmetric, have been used to investigate the model in [17, 18]. Since the Thirringmodel is not a gauge theory, the long-standing objections [19] to the SLAC approachassociated with lack of localisation of the fermion-gauge vertex do not apply. No evidencehas been found for spontaneous symmetry breaking for any integer value of N (ie. anyunitary local version of the model); an estimate N c = 0 . < N ) symmetry is notmanifest 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. Anaspect worth highlighting is that in the most promising approach the fermions interactwith the auxiliary field A µ throughout the bulk, ie for 0 < x < L s , very similar to theway gauge theories are formulated with DWF [21].We have found that the L s → ∞ limit becomes particularly challenging preciselyin the strong-coupling regime where symmetry breaking is anticipated. In particular,Ref. [14] presented results for N = 0 , , × L s (for N >
0) with L s ranging from8 to 40. The results for N = 0 and N = 2 are fairly clear-cut: symmetry-breakingis observed in the former case and not in the latter. For N = 1 there are qualitativeindications of a change in the system’s behaviour at the strongest coupling explored( g − a = 0 .
3, where the lattice spacing a defines the scale). In particular the bilinearcondensate (cid:104) ¯ ψψ ( m ) (cid:105) displays significant L s -dependence in this regime, and the resultingsignal estimated in the L s → ∞ limit is significantly greater than that from weakercouplings; there is also a marked departure from linear m -dependence. The conclusionreached in [14] is that 0 < N c < N c > N = 1 is numerically more demanding [14]. Another is that in afinite volume the bilinear condensate vanishes identically for m = 0 so that neither orderparameter nor associated susceptibility are directly accessible in the U(2 N )-symmetriclimit. In the current paper we apply improved code and substantially enhanced com-puting resources to both issues, by exploring the model on both 12 × L s and 16 × L s ,with L s as large as 48, with much finer resolution along the coupling axis, particularlyin the suspected critical region. The resulting order parameter data (cid:104) ¯ ψψ ( g , m ) (cid:105) areobtained with sufficient statistical precision to permit fits to a renormalisation-group in-spired equation of state (EoS) applicable away from m = 0, yielding estimates for bothcritical coupling g 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 > β ≡ g − a , m , L s and space-time volume,in the regime 0 . ≤ β ≤ .
52. Results for (cid:104) ¯ ψψ (cid:105) are first extrapolated to L s → ∞ ,and then fitted to the empirical EoS (16) below. The picture that emerges appearsunder control; the L s -extrapolation and EoS-fitting procedures almost, but not quite,commute, yielding fairly stable estimates for β c ≈ .
28 and estimates for the exponents β m , δ distinct from those of mean field theory. Finite volume effects appear remarkablysmall. However, all the data used in this analysis lie in the symmetric phase. To addressthis Sec. 3.2 presents a complementary study at fixed L s = 48, but this time exploring4ouplings as strong as β = 0 .
23. Both (cid:104) ¯ ψψ (cid:105) and its associated susceptibility χ (cid:96) arestudied. Here some issues emerge; the EoS fit is not so successful at describing datain the strong-coupling phase, and χ (cid:96) , 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(2 N ) symmetry, and the bose action den-sity (2 g ) − (cid:104) A µ (cid:105) . In Sec. 3.3 we present results for the locality, and recovery of theGinsparg-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 localU(2 N )-symmetric field theory at the QCP. Our conclusions are discussed in Sec. 4. We use an auxiliary field formulation to represent the Thirring interaction, which re-casts the fermion action as a bilinear form while preserving the global symmetries. Incontinuum notation the Lagrangian density reads L (cid:48) = ¯ ψ ( ∂/ + iA/ + m ) ψ + 12 g A µ . (6)On a lattice the vector auxiliary field A µ is naturally formulated on a link. Preservingfermion global symmetries while transcribing to a lattice is of course a long-standingissue in lattice field theory. Our approach uses domain wall fermions (DWF), in whicha fictitious third dimension is introduced so the fermions Ψ( x, s ) are defined in 2+1+1 d : S = (cid:88) x,y (cid:88) s,s (cid:48) ¯Ψ( x, s )[ δ s,s (cid:48) D W x,y [ A ] + δ x,y D s,s (cid:48) ]Ψ( y, s (cid:48) ) + mS m + 12 g (cid:88) x,µ A µ ( x ) . (7)We use the M¨obius formulation, implying that D W [ A = 0] , D are free Wilson derivativeoperators, with D s,s (cid:48) having open boundary conditions implemented on the hoppingterms, viz. δ s ∓ ,s (cid:48) (1 − δ s (cid:48) , /L s ), 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 modelemploys 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 generatessolely 4-fermi terms, and not higher point contact interactions. Further details are setout in [14].In the large- L s limit the free kinetic operator approaches an overlap operator D OL defined in 2+1 d and satisfying a generalisation of the Ginsparg-Wilson relations [23, 24] { γ , D OL } = 2 D OL γ D OL ; { γ , D OL } = 2 D OL γ D OL ; [ γ γ , D OL ] = 0 . (8)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 ofthe continuum model. This only holds for the bulk formulation of the Thirring model5ith DWF, and not the alternative “surface” formulation discussed in [20, 14]. Forstrongly-interacting fields the recovery of the GW relations (8) and the locality of thecorresponding D OL will be discussed further in Sec. 3.3 below.The bridge between 2+1+1 d and the target model rests on the identification ofphysical fermion fields in 2+1 d localised entirely on the domain walls, which we regardas a working assumption: ψ ( x ) = P − Ψ( x,
1) + P + Ψ( x, L s ); ¯ ψ ( x ) = ¯Ψ( x, L s ) P − + ¯Ψ( x, P + , (9)with projectors P ± ≡ (1 ± γ ). 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 ¯ ψγ ψ and im ¯ ψγ ψ . For DWF fermions with L s finite the choice im ¯ ψγ ψ with ψ, ¯ ψ specified in (9) gives the best approach to L s → ∞ [22, 24, 14], and we usethis definition throughout this paper. Approach to the U(2) symmetric limit will bemonitored in Sec. 3.2 below via the residual δ h (cid:39) (cid:104) ¯ ψψ (cid:105) − (cid:104) ¯ ψiγ ψ (cid:105) [22].Specifying the bilinear component of the action (7) by M , then the action simulatedusing bosonic pseudofermion fields Φ , Φ † is S pf = Φ † (cid:110) [ M † M m h =1 ] [ M † M m = m ] − [ M † M m h =1 ] (cid:111) Φ . (10)Assuming det M is real, then the resulting functional weight det[ M m = m M − m h =1 ] tendsto det D OL ( m ) in the limit L s → ∞ [24]. An RHMC algorithm is needed to handlethe fractional powers in (10), as described in [14]. Some tests of the accuracy of therational approximation needed to implement fractional powers in the parameter regimeof interest are presented in Sec. 3.2. In the meantime the code has been modified intwo main aspects: multiprocess parallelisation and a simplified, more stable solver. Theparallelisation makes use of the MPI paradigm, with a 3D domain decomposition acrossthe spatial and temporal directions, leaving the Domain Wall dimension uncut to allowthe 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 QMR-based one, and being stable in single precision. The possibility of running the solver insingle precision during the molecular dynamics part of the algorithm leads to another50% save in memory which, as the solver is severely memory-bound, translates into adirect performance increase. We have made the simulation code publicly available [25].Finally, the observables studied in the bulk of the paper are simply the bilinearcondensate (cid:104) ¯ ψψ (cid:105) ≡ ∂ ln Z ∂m ≡ (cid:104) Σ (cid:105) (11)and its associated susceptibility χ (cid:96) = (cid:104) Σ (cid:105) − (cid:104) Σ (cid:105) . (12)The notation is slightly loose; the bilinear actually used in computations should beunderstood to be i ¯ ψγ ψ . The susceptibility χ (cid:96) defined in (12) is often referred to as the6 isconnected component, and is expected to manifest any singular behaviour expectedat a continuous phase transition. The full susceptibility corresponding to ∂ ln Z /∂m contains an extra connected component not included in the variance of the bilinear orderparameter calculated here. In our work the expectations (11,12) are calculated withunbiassed estimators using 10 independent Gaussian noise vectors per configuration. , , and the Extrapolation L s → ∞ In previous work [14], simulations of the N = 1 model on a 12 × L s system found signalsconsistent with broken U(2) symmetry in the vicinity of β ≡ g − a ≈ .
3. Crucially, inorder to demonstrate this it proved necessary to explore a large range 8 ≤ L s ≤ L s → ∞ . In this subsection we revisit this issue,this time exploring the coupling range 0 . ≤ β ≤ .
52 with much finer resolution,with fermion masses ma = 0 . , . . . , .
05, spacetime volumes 12 and 16 and L s =8 , , . . . ,
48. The data presented result from at least 3000 RHMC trajectories of meanlength 1.0 for each parameter set { β, m, L s } .First let’s discuss the L s → ∞ extrapolation; empirically the convergence to thelarge- L s limit is exponential [14]: (cid:104) ¯ ψψ (cid:105) L s = ∞ − (cid:104) ¯ ψψ (cid:105) L s = A ( β, m ) e − ∆( β,m ) L s . (13)Sample fits are shown in Fig. 1. The extracted decay constant ∆( β, m ) is shown inFig. 2. While errors are appreciable, particularly at weaker couplings where the signalis small, it can be seen that for the weaker couplings 0 . ≤ β ≤ .
44 ∆ ≈ .
03 and isapproximately m -independent, whereas for the stronger couplings 0 . ≤ β ≤ .
34, ∆( m )is roughly linear, with an intercept ∆( m = 0) ≈ . L s ∼ O (10 )is needed to completely control the extrapolation in this regime.Next we turn to analysis of the critical point, assumed present at ( β, m ) = ( β c , m (cid:54) = 0, wefollow an indirect route, motivated by the renormalisation group [7]. At fixed spacetimevolume, assume a scaling form m ( (cid:104) ¯ ψψ (cid:105) , t ) = (cid:104) ¯ ψψ (cid:105) δ F ( t (cid:104) ¯ ψψ (cid:105) − βm ) (14)where t ≡ β − β c . In the limits m = 0, t = 0 we immediately recover | t | β m ∝ (cid:104) ¯ ψψ (cid:105) ; m = F (0) (cid:104) ¯ ψψ (cid:105) δ (15)whereupon we identify the conventional critical exponents β m and δ familiar from theferromagnetic transition. For small t we may Taylor-expand the scaling function F toyield the equation of state (EoS): m = A ( β − β c ) (cid:104) ¯ ψψ (cid:105) δ − βm + B (cid:104) ¯ ψψ (cid:105) δ . (16)7
10 20 30 40 50 L s . . . . . . . ¯ ψψ β : 0 . β : 0 . β : 0 . β : 0 . β : 0 . β : 0 . β : 0 . β : 0 . β : 0 . β : 0 . β : 0 . β : 0 . β : 0 . Figure 1: Bilinear condensate (cid:104) ¯ ψψ ( L s ) (cid:105) for ma = 0 .
01 on 16 together with a fit to (13). .
00 0 .
01 0 .
02 0 .
03 0 .
04 0 . am . . . . . . . . ∆ β : 0 . β : 0 . β : 0 . β : 0 . β : 0 . β : 0 . β : 0 . β : 0 . Figure 2: The decay constant ∆( β, m ) defined in (13) on 16 , showing distinct trends in the couplingranges β ≥ .
36 and β ≤ .
34. Data points have been shifted horizontally for improved readability. .20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.600.000.050.100.150.200.25 /dof=27.2 m Ndat = 72 m=0m=0.005m=0.01m=0.02m=0.03m=0.04m=0.05
Figure 3: Fit to the EoS (16) for 16 ×
48. The full line shows the fit extrapolated to the U(2)-symmetriclimit m = 0. The large χ value is dominated by data with β = 0 .
30 – see Fig. 7 below.
Results of a least-squares fit of (16) to data from 16 ×
48, 0 . ≤ β ≤ .
52, and ma = 0 . . . . , .
05 are shown in Fig. 3, together with the curve resulting fromthe same fit parameters as m →
0. The data support a symmetry-breaking continuousphase transition at β c = 0 . χ / dof results from excluding the data with β = 0 .
3; here we retain them for the sake ofuniformity in the subsequent analysis). The fitted exponents β m = 0 . δ = 3 . β m = , δ = 3.However, in order to identify this transition with the breaking of U(2) symmetrya minimum requirement is stability under the extrapolation L s → ∞ . The purist’sapproach, 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 onthe datapoints (and correspondingly smaller χ / dof) – the extrapolation is particularlydifficult to control at weak coupling where the signal is small. The fitted EoS stillsupports a continuous phase transition, though with modified exponents β m = 0 . δ = 4 . β c = 0 . ma = 0 . , .
01; however, due to the large errorbars exclusion of thesemasses does not significantly alter the fit.A more pragmatic approach is to perform fits at fixed L s and then study the be-9 .20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.600.000.050.100.150.200.25 /dof=2.22 m Ndat = 72
Figure 4: Fit to the EoS (16) for 16 extrapolated to infinite L s . haviour of the fit parameters as L s → ∞ . This is less well-motivated theoretically, buthas the practical advantage that the data passed to the least-squares fitting procedureis of much higher statistical quality. Results for the critical coupling and the exponentsare shown in Figs. 5, 6 respectively. In Fig. 5 it is notable that β c ( L s ) from 12 and 16 volumes are compatible; the evolution with L s is smooth but falls significantly short ofthe fit of the extrapolated data even by L s = 48, so it is not clear whether extrapolationand fitting commute. In Fig. 6 by contrast the fitted β m ( L s ) and δ ( L s ), while agreeingat large L s , approach this limit from opposite directions on 12 and 16 ; moreover whileon 12 the data from the largest available L s are compatible with the values extractedfrom the extrapolated data, as already remarked there is a significant disparity on 16 . , L s = 48 The results of the previous subsection support the claim made in [14] that the N = 1model has a continuous phase transition associated with the onset of a bilinear fermioncondensate for β ≈ .
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 withcritical exponents both distinct from mean-field values, and consistent with a QCPdefining a previously unknown strongly coupled quantum field theory of fermions. Wemight be concerned, however, that all data analysed in Sec. 3.1 lie on the weak-couplingside of the putative transition, and that there are hints from Fig. 3 and particularly10
10 20 30 40 50 L s β c Figure 5: Fitted critical coupling β c on both 12 and 16 spacetime volumes as a function of L s . Thefull (16 ) and dashed (12 ) lines show the values from fits of L s → ∞ extrapolated data. β m MFT L s δ Figure 6: Fitted critical exponents β m (upper panel) and δ (lower panel) on both 12 and 16 spacetimevolumes as a function of L s . Dash-dotted lines show mean field values .20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.600.000.050.100.150.200.25 /dof=5.67 m Ndat = 66 m=0m=0.005m=0.01m=0.02m=0.03m=0.04m=0.05
Figure 7: The fixed- L s dataset used for the analysis of Sec. 3.2, together with a fit to the EoS (16).The fit parameters are β c = 0 . β m = 0 . δ = 3 . Fig. 4 that the EoS fit is not doing such a good job as either m → β → β c + .To investigate further we now turn to data generated using approximately 5000 RHMCtrajectories at fixed L s = 48 on spacetime volume 16 , but extending to β -values on thestrong-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 (cid:104) ¯ ψψ ( β, m ) (cid:105) dataset, and a fit to (16) based on 0 . ≤ β ≤ . β = 0 .
30. It is immediately apparent, however, that itfails to model data on the strong coupling side of the transition; here (cid:104) ¯ ψψ (cid:105) falls belowthe model, the effect is more pronounced with decreasing m , until by ma = 0 .
005 thecurve becomes flat for β < ∼ .
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 β → + . A possible origin of this behaviour was suggested in [7]: as a result of the linearcoupling between A µ and the fermion current, the leading order large- N correction to theauxiliary propagator is not transverse in a lattice regularisation, leading to breakdownof positivity as β →
0. In the large- N approach the effect is mitigated by an additiverenormalisation of g − , so that the strong coupling limit of the model is now taken as β → β ∗ >
0. In [11] β ∗ was taken to be the location of the maximum of (cid:104) ¯ ψψ ( β ) (cid:105) ,enabling a model equation of state (cid:104) ¯ ψψ ( N ) (cid:105) in the effective strong-coupling limit and12 .20 0.25 0.30 0.35 0.40 0.45 0.50 0.550.000.010.020.030.040.050.06 N pf = 25 N pf = 29 L s = 64 Figure 8: (cid:104) ¯ ψψ ( β ) (cid:105) for ma = 0 .
005 on 16 ×
48. The dashed line is the same EoS fit shown in Fig. 7.Also shown is the result of a pilot simulation with β = 0 . L s = 64. consequent prediction of N c . A decrease of (cid:104) ¯ ψψ (cid:105) β → has also been reported in simulationswith SLAC fermions [17, 18] and with DWF in a variant “surface” formulation of themodel [20], suggesting that strong coupling lattice artifacts are a generic feature of theThirring 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 separationbetween β ∗ and any β c associated with a Thirring model QCP. From Fig. 8 we mightestimate β ∗ ≈ .
25, uncomfortably close, with current resolution, to the β c estimates ofSec. 3.1.At this point it is appropriate to discuss a technical aside. In the RHMC algorithmdescribed 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 A p (cid:39) r p ( A ) = α + N pf (cid:88) i =1 α i A + β i , (17)where the coefficients α i , β i may be calculated using the Remez algorithm implemen-tation described in [26]. They are chosen so that over a spectral range ( λ d , . | r p ( x ) − x p | < − for matrices needed during trajectory guidance and < − for thoseneeded in the Monte Carlo acceptance step. For all work to date we have used λ d = 10 − corresponding to the smallest value of ( ma ) explored, which translates to partial fractionnumbers N pf = 12 (guidance) and N pf = 25 (acceptance); however one might question13hether this is sufficiently accurate for studies with ma = 0 . β -values with Remez coefficients generatedwith λ d = 10 − , 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 datacalculated using the previous λ d = 10 − .Next we present data for the susceptibility χ (cid:96) defined in (12), for the whole datasetin Fig. 9 and for the lightest ma = 0 .
005 in Fig. 10. As might be anticipated, statisticalerrors in χ (cid:96) are considerably larger than those for the condensate, and accordingly wechoose not to attempt an L s → ∞ extrapolation. However, again, the agreement be-tween results obtained using the default and enhanced rational approximations seen inFig. 10 is reassuring. For each value of m χ (cid:96) ( β ) is non-monotonic, the peak shifting tostronger coupling as m decreases in accord with expectations for a second derivative ofthe free energy at a critical point; this is corroborated by the model prediction obtainedby differentiation of (16) with respect to m , and plotted using the fitted parameters inFig. 11. Fig. 10 suggests that the location and even the sharpness of the peak at criti-cality is roughly as expected, once an empirical rescaling is applied. However, there arefeatures of Fig. 9 which are clearly problematic; the m -ordering of the data is oppositeto model expectations, with χ (cid:96) increasing with m over the whole β -range studied, andthe convergence of χ (cid:96) 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 asfollows: δ h ( L s ) = Im (cid:104) ¯Ψ( x, s = 1) iγ Ψ( x, s = L s ) (cid:105) = − Im (cid:104) ¯Ψ( x, s = L s ) iγ Ψ( x, s = 1) (cid:105) . (18)In [22] 2 δ h was found to furnish a lower bound for the difference between (cid:104) ¯ ψψ (cid:105) and (cid:104) ¯ ψiγ ψ (cid:105) , and to vanish ∝ e − cL s for quenched QED . In the Thirring model, δ h ( L s ) forvarious couplings was presented in Fig. 12 of [14]. While in all cases δ h still decreaseswith L s , it grows larger as coupling increases, and by β = 0 . c evendevelops a dependence on m . Fig. 12 taken at fixed L s = 48 confirms that m -dependenceof δ h does indeed set in for β < ∼ .
4, and that δ h continues to grow as β decreases,suggesting that recovery of U(2) symmetry will be an ever-increasing challenge in thesymmetry broken phase as m → g ) − (cid:104) A µ (cid:105) .As discussed in [14], for DWF with the bulk formulation of the Thirring model thereis no simple interpretation in terms of a local four-fermion condensate available; ratherwe regard it as an extra observable sensitive to light fermion dynamics. Its behaviour isnon-monotonic, with a minimum at β (cid:39) .
46 before rising to approach and then exceedthe free-field value at β (cid:39) .
24. The notable feature of Fig. 13 is the fermion mass-dependence; broadly speaking the departure from the free-field result increases withdecreasing m (although the m -ordering of the data is somewhat noisy), the effect beingmost pronounced for 0 . < ∼ β < ∼ . .20 0.25 0.30 0.35 0.40 0.45 0.50 0.550.00.20.40.60.8 m=0.005m=0.01m=0.02m=0.03m=0.04m=0.05 Figure 9: Susceptibility χ (cid:96) ( β, m ) on 16 × N pf = 25 N pf = 29 L s = 64 Figure 10: Susceptibility χ (cid:96) ( β ) for ma = 0 .
005 on 16 ×
48. The dashed line is calculated using thesame EoS fit shown in Fig. 7, multiplied by an empirical factor 0.014. Also shown is the result of a pilotsimulation with β = 0 . L s = 64. .2 0.3 0.4 0.5 b c l ma= ma= ma= ma= ma= ma= Figure 11: Susceptibility χ (cid:96) ( β, m ) obtained by differentiation of the EoS (16). The inset shows the“corrected” version discussed in Sec. 4 h m=0.005m=0.01m=0.02m=0.03m=0.04m=0.05 Figure 12: The U(2)-breaking residual δ h ( β, m ) on 16 × .20 0.25 0.30 0.35 0.40 0.45 0.50 0.551.3251.3501.3751.4001.4251.4501.4751.5001.525 A / g free-field result m=0.005m=0.01m=0.02m=0.03m=0.04m=0.05 Figure 13: Auxiliary action density (2 g ) − (cid:104) A µ (cid:105) vs. β on 16 ×
48. The dashed line through the ma = 0 .
005 data is merely to guide the eye.
The equivalence of DWF [27, 21] and the (truncated) overlap operator [28] is well es-tablished in 3+1 d , eg. [29]. This equivalence is further shown in 2+1 d [24] for both theregular mass term m ¯ ψψ and the linearly independent twisted mass terms im ¯ ψγ , ψ intro-duced above. As such, locality of the domain wall operator in the target dimensionalitycan be demonstrated by showing the locality of the overlap operator.We use the Shamir and Wilson formulations of the overlap operator with twistedmass − im ¯ ψγ ψ given by D OL ( m ) = 1 − imγ imγ V S/W (19)where the Wilson and Shamir kernels, defined via V = γ sgn( H ), are H W = γ D W H S = γ D W D W (20)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 relation17 γ = ( γ H ) † is preserved. Our formulation of DWF in the L s → ∞ limit is expectedto recover (19) with the Shamir kernel [24].The standard mass formulation of the overlap operator is D IOL ( m ) = m + − m V S/W .The signum function typically is approximated with a rational function resulting in atruncated overlap operator. The hyperbolic tangent (polar) approximation is the mostcommonly used and may be expressed as [30]sgn( x ) ≈ tanh( n tanh − x ) = xn (cid:81) n/ − j =1 [ x + (tan jπn ) ] (cid:81) n/ − j =0 [ x + (tan ( j +1 / πn ) ] (21)for n even. For the Shamir formulation it is much more efficient to use a formulationexploiting an extra dimension, and the truncated overlap operator can be reconstructeddirectly from a domain wall formulation. The standard mass and alternative massdomain wall operators may be expressed as D I/ DW = D DW + mD I/ DW , where (for n ≡ L s = 4) D DW = D W + I − P − − P + D W + I − P − − P + D W + I − P − − P + D W + I (22) D IDW = P + P − , D DW = iγ P + iγ P − (23)Then with C = P − P + P − P +
00 0 P − P + P + P − , C † = C − = P − P + P + P − P + P −
00 0 P + P − (24)we have the following relation [24, 29], where the precise form of the (cid:52) terms is unim-portant for our purposes. K = C † ( D IDW (1)) − D DW ( m ) C = D OL ( m ) 0 0 0 − (1 − m ) (cid:52) R − (1 − m ) (cid:52) R − (1 − m ) (cid:52) R (25)So the overlap operator (19), with the Shamir kernel (20), truncated with the polarapproximation (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 IDW (1)) − D IDW ( m ) C .18nvariant transformations, corresponding to the Ginsparg-Wilson (GW) relations (8)are given by [22] Ψ → e iαγ (1 − aD ) Ψ ; ¯Ψ → ¯Ψ e iαγ (1 − aD ) Ψ → e iαγ (1 − aD ) Ψ ; ¯Ψ → ¯Ψ e iαγ (1 − aD ) Ψ → e iαγ γ Ψ ; ¯Ψ → ¯Ψ e iαγ γ (26)This relation is exact for the 2+1 d overlap operator D = D OL , and is reproduced byDWF in 2+1+1 d in the L s → ∞ limit. In order to recover the U (2) symmetry in thecontinuum limit a →
0, we must have the GW terms aDγ , D and equivalently thetransform terms aD vanishing in the same limit. A sufficient condition for this to be thecase is the Dirac operator being exponentially local. Evidence for this has been givenfor the overlap operator in 3+1 d [31], and it behooves us to investigate the 2+1 d case.Since the overlap operator is a dense matrix and manifestly non-local, demonstration ofexponential 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 e iαγ (1 − aD ) Ψ = ( I + iαγ (1 − aD · · · )Ψ (27)so that recovery requires [ aD Ψ] a → = 0 (28)We have Ψ (cid:48) j = [ a (cid:80) i D ji Ψ i ] a → = 0 which is true if [ (cid:80) i D ji Ψ i ] a → < ∞ which is true forany bounded Ψ if [ (cid:80) i D ji ] a → < ∞ which is true if D is exponentially local, and henceexponential locality allows recovery of U (2) symmetry.In order to illustrate the locality of the overlap operator, we follow [31]. Let ψ ( x ) = Dη ( x ) (29)where the point source η ( x ) = δ x,y δ α, where y is an arbitrary location and α a spinorindex. Then we calculate the decay as f ( r ) = max {|| ψ ( x ) || : || x − y || = r } (30)where the “Manhattan taxi distance”, || x − y || = (cid:80) µ | x µ − y µ | is just the L norm.The locality of D OL in the critical region is illustrated for the twisted im ¯ ψγ ψ massterm with ma = 0 .
005 in Figs. 14,15. Within the limitations imposed by a 16 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 compar-ison are data obtained with unitary link fields e ± iA µ , where exponential localisation ismore manifest. Convergence to a meaningful value of the decay rate on this small volumeis seen to be difficult in Fig. 15 showing the decrement from one r -value to the next, byplotting f ( r ) /f ( r − r f(r) β =0.24 β =0.27 β =0.30 β =0.33 β =0.30 unitary Figure 14: Localisation of the overlap D with Shamir kernel and L s = 48 averaged over 32 sources oneach of 5 configurations r f(r)/f(r-1) β =0.24 β =0.27 β =0.30 β =0.33 β =0.31 unitary Figure 15: The decrement f ( r ) /f ( r −
1) for the data of Fig. 14 L s δ GW β =0.30 S β =0.33 S β =0.30 W β =0.33 W Figure 16: GW error δ GW calculated on a single configuration for Shamir (S) and Wilson (W) kernels.The inset compares results obtained using unitary link fields at β = 0 . marginal variation, indicating locality is not jeopardised at a critical point. Varying m does not change the conclusions, nor does using the Wilson kernel rather than Shamir.We also examine the truncation/finite- L s error of the overlap/DWF operator via theGW term, as a means to assess recovery of U(2) symmetry. In the n → ∞ ( L s → ∞ )limits the GW error δ GW , given as δ GW = || ( γ D + Dγ − Dγ D ) φ || ∞ (31)with φ a complex field chosen at random at each lattice site such that each componentis distributed uniformly in ( − , δ GW in Fig. 16. Although the boson field is generatedwith 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 similarto δ h defined in (18) (plotted as a function of L s in Fig. 12 of [14]). It is interestingthat convergence with the Wilson kernel is slightly faster as compared to the Shamirformulation, which is surprising since with non-unitary links the Wilson kernel is notbounded which should prejudice convergence [30, 24]. Again, for comparison resultsobtained with unitary link fields are included in the inset; here exponential falloff of theerror is much sharper. These results suggest that the very large values of L s needed forU(2 N ) recovery in the critical region have their origin in the non-unitary nature of thelink fields in this formalism. 21 Discussion
The main conclusion of this study is that we can state with much more confidence thatthere is a phase transition associated with bilinear condensation in the Thirring modeldefined with N = 1 domain wall fermions, and hence that N c > L s → ∞ extrapolation, enabling fits to an RG-inspired equation of statewith β c (cid:39) . β > ∼ β c . The fitted critical exponentsare significantly different from their mean-field values, although both Fig. 2 and Figs. 5,6suggest simulations with larger L s , perhaps even L s ∼ O (10 ), will be needed in orderto fully control the predictions. Further encouraging signs are a susceptibility peak ofthe expected shape seen in Fig. 10 and the m -dependence of the bosonic action in thesubcritical region seen in Fig. 13.The result N c > ∼ N c =0 . A µ resembles a gauge field. Ultimately, physics at a QCP is specified not by a Lagrangiandensity, either continuum-like or regularised, but by more primitive considerations suchas dimensionality, field content, and of course the pattern of global symmetry breaking.Interestingly, a recent study using the Conformal Bootstrap predicts 1 < N c < [32].For the first time in Sec. 3.3 we have studied properties of the 2 + 1 d overlap operator D OL associated with DWF in the vicinity of the critical point. Fig. 14 suggests thereexists a limit where aD OL →
0, a necessary condition for the existence of a continuumlimit with conventional U(2 N ) symmetry, although confirmation will require studies onlarger spacetime volumes. However, symmetry recovery also requires the restoration ofthe Ginsparg-Wilson relations (8). Both the direct estimates shown in Fig. 16 and theindirect measure via the residual δ h shown in Fig. 12 suggest this is at best restoredonly very slowly in the critical region. Note that δ h → O ( e − ∆ L s ) corrections are applied to the order parameter as described in Sec. 3.1. Itis now becoming apparent that this behaviour has its origins in the use of non-unitarylink fields in the Wilson operator D W ; recall this originates in the desire to have onlyfour-point fermion interactions left once the auxiliary is integrated over. Fig. 16 is alsoa motivation to explore measurement using D OL defined with the alternative Wilsonkernel, and/or with expansion order n greater than the L s used in ensemble generationwith 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 theexpectation of Fig. 11, and moreover requires a seemingly arbitrary rescaling even to fit22 single m value. The behaviour at weaker couplings can be somewhat accommodatedby replacing the LHS of the EoS (16) by a factor m ( m/m ) α . The inset of Fig. 11 showsthe susceptibility curves thus obtained with m a = 0 . α = 0 .
3. This fails to fix thehierarchy in the critical region, however; it seems safer to conclude that χ (cid:96) defined in(12) is not the second derivative of the free energy of a 2+1 d theory. A modificationof the physical field prescription (9) may be needed, to allow for the possibility thatthe relevant modes close to the walls “leak” somewhat into the bulk as ma →
0. Arelated puzzle, again a failure to reconcile the observed behaviour of order parameterand susceptibility data with possibly the same cause, is the breakdown of the axial WardIdentity noted in [20].The flattening of the order parameter at small β seen in Figs. 7,8 remains a worryingaspect. 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 ×
64 with β = 0 . ma = 0 . L s = 64 point lying on the fitted curve in Fig. 8should be regarded at this stage as a coincidence). Another possibility is that the scalingwindow where (16) is applicable is simply very narrow on the horizontal scale used inthese figures. As things stand, however, we have not yet demonstrated a clear separationbetween β c and a putative β ∗ where lattice artifacts dominate, and hence do not yethave an understanding of the broken phase comparable with that for β > β c . The slowconvergence to the U(2 N ) limit discussed above may prove a formidable obstacle; thereis much work still to be done. Acknowledgements
Sunbird facility of Supercom-puting Wales. The work of SJH was supported by STFC grant ST/L000369/1, of MMin part by the Supercomputing Wales project, which is part-funded by the EuropeanRegional Development Fund (ERDF) via Welsh Government, and of JW by an EPSRCstudentship.
References [1] M. Gomes, R.S. Mendes, R.F. Ribeiro and A.J. da Silva, Phys. Rev. D (1991)3516.[2] T. Itoh, Y. Kim, M. Sugiura and K. Yamawaki, Prog. Theor. Phys. , 417 (1995).233] M. Sugiura, Prog. Theor. Phys. (1997) 311.[4] H. Gies and L. Janssen, Phys. Rev. D (2010) 085018;[5] L. Janssen and H. Gies, Phys. Rev. D (2012) 105007;[6] F. Gehring, H. Gies and L. Janssen, Phys. Rev. D (2015) 085046.[7] L. Del Debbio, S. Hands and J.C. Mehegan, Nucl. Phys. B (1997) 269.[8] S. Kim and Y. Kim, hep-lat/9605021. [9] L. Del Debbio and S.J. Hands, Nucl. Phys. B , 339 (1999).[10] S. Hands and B. Lucini, Phys. Lett. B (1999) 263.[11] S. Christofi, S. Hands and C. Strouthos, Phys. Rev. D (2007) 101701(R).[12] S. Chandrasekharan and A. Li, Phys. Rev. Lett. (2012) 140404.[13] C. Burden and A.N. Burkitt, Europhys. Lett. (1987) 545.[14] S. Hands, Phys. Rev. D (2019) no.3, 034504[15] S. Hands, Phys. Rev. D (1995) 5816.[16] T. Appelquist, A.G. Cohen and M. Schmaltz, Phys. Rev. D (1999), 045003.[17] B.H. Wellegehausen, D. Schmidt and A. Wipf, Phys. Rev. D (2017) 094504[18] J.J. Lenz, B.H. Wellegehausen and A. Wipf, Phys. Rev. D (2019) no.5, 054501.[19] L.H. Karsten and J. Smit, Phys. Lett. B (1979), 100-102.[20] S. Hands, JHEP (2016) 015.[21] Y. Shamir, Nucl. Phys. B (1993), 90-106.[22] S. Hands, JHEP (2015) 047.[23] P.H. Ginsparg and K.G. Wilson, Phys. Rev. D (2016) 264.[25] E. Bennett, S. Hands, M. Kappas and M. Mesiti, https://doi.org/10.5281/zenodo.4016827 [26] M.A. Clark and A.D. Kennedy, https://github.com/mikeaclark/AlgRemez ,(2005).[27] D.B. Kaplan, Phys. Lett. B (1992), 342-347.2428] H. Neuberger, Phys. Lett. B (1998), 141-144;Phys. Lett. B (1998), 353-355.[29] R.C. Brower, H. Neff and K. Orginos, Comput. Phys. Commun. (2017), 1-19.[30] A.D. Kennedy, Algorithms for dynamical fermions , hep-lat/0607038 .[31] P. Hernandez, K. Jansen and M. L¨uscher, Nucl. Phys. B (1999), 363-378[32] Z. Li, arXiv:1812.09281 [hep-th]arXiv:1812.09281 [hep-th]