The Phase diagram of the three dimensional Thirring model

arXiv:hep-lat/9906008v1 11 Jun 1999 SWAT/99/231 IFUP-TH 31/99 June 1999 The Phase Diagram of the Three Dimensional Thirring Model Simon Hands a a and Biagio Lucini b,c Department of Physics, University of Wales Swansea, Singleton Park, Swansea SA2 8PP, U.K. b c Scuola Normale Superiore, Piazza dei Cavalieri 7, I-56126 Pisa, Italy. INFN Sezione di Pisa, Via Vecchia Livornese 1291, I-56010 S. Piero a Grado (Pi), Italy. Abstract We present Monte Carlo simulation results for the three dimensional Thirring model on moderate sized lattices using a hybrid molecular dynamics algorithm which permits an odd or non-integer number Nf of fermion flavors. We find a continuous chiral symmetry breaking transition for Nf ≃ 3 with critical exponents consistent with expectations from previous studies. For Nf = 5 the order of the transition is difficult to determine on the lattice sizes explored. We present a phase diagram for the model in the (1/g 2 , Nf ) plane and contrast our findings with expectations based on approximate solutions of the continuum Schwinger-Dyson equations. PACS: 11.10.Kk, 11.30.Rd, 11.15.Ha Keywords: four-fermi, Monte Carlo simulation, dynamical fermions, chiral symmetry breaking, renormalisation group fixed point 1 1 Introduction and Algorithm Recent numerical investigations [1, 2, 3] have revealed an interesting phase structure for the three dimensional Thirring model as a function of coupling g 2 and number of fermion species Nf . For sufficiently small Nf < Nf c , there is a continuous transition at g 2 = gc2 (Nf ) between a weak coupling phase in which chiral symmetry is realised in the limit m → 0, and a strong coupling one in which the symmetry is spontaneously broken. The critical indices characterising the transition are distinct for Nf = 2 and Nf = 4, suggesting that the continuum limits defined at the critical points define distinct interacting field theories. For Nf = 6, however, simulations on a 163 lattice with m = 0.01 give tentative evidence for a first order chiral transition [3], implying both that no continuum limit exists in this case, and that 4 < Nf c < 6. Both the nature of the transition and the value of Nf c are non-perturbative issues, inaccessible via either a standard perturbative expansion (which is non-renormalisable for d > 2), or a 1/Nf expansion [4, 5]. There have, however, been analytic attempts to investigate the transition via the truncated Schwinger-Dyson (SD) equations [4, 6, 7, 8]. In the most systematic treatment [6], the SD equations are solved in ladder approximation in the the strong coupling limit; chiral symmetry breaking solutions are found for Nf < Nf c = 128/3π 2 ≃ 4.32. Note that in a continuum approach the value of Nf c manifests itself as the value of Nf below which a non-trivial solution can be found, eg. by bifurcation theory. In a later paper [8], Sugiura extended the solution to finite g 2 for d ∈ (2, 4); for the case d = 3 his predictions for the critical coupling as a function of Nf read: g2 ≪ 1 : gc2 = g2 ≫ 1 : gc2 = with 2π 2 Nf 3  (1)  1 2π − 4 tan−1 ω − 1 ω Nf c −1 ω 2 (Nf ) = Nf  6 exp  (2) The essentially singular behaviour seen in (2) as Nf ր Nf c is consistent with the existence of a conformal fixed point [9]. Using a different sequence of truncations of the SD equations, however, Hong and Park [7] found chiral symmetry breaking solutions for all Nf , and predicted gc2 ! π2 ∝ exp Nf , 16 2 (3) which is manifestly non-perturbative in 1/Nf . The differing predictions of the SD approach motivate the use of lattice field theory methods to address this problem. In this Letter we present results of simulations performed with values of Nf falling between the values Nf = 2, 4 and 6 explored in previous works, in an effort both to map out the phase diagram in more detail, and to constrain further the value of Nf c . To start with, let us define the lattice action: 1X S = χ̄i (x)ηµ (x) (1 + iAµ (x))χi (x + µ̂) − (1 − iAµ (x − µ̂))χi (x − µ̂) 2 xµi  N X 2 A (x) 4g 2 xµ µ xi X N X 2 χ̄i (x)Mij [A](x, y)χj (y) + 2 ≡ A (x), 4g xµ µ xyi +m X χ̄i (x)χi (x) +  (4) where the indices i, j run over N flavors of staggered lattice fermion. The Aµ are real auxiliary vector fields defined on the lattice links; gaussian integration over Aµ yields a form of the action with explicit four-fermion couplings [2, 3] (other variants of the lattice-regularised Thirring model are discussed in [2]). Analysis of the spin-flavor content of staggered lattice fermions [10] reveals that in 3 dimensions the number of continuum four-component physical flavors is given by Nf = 2N. Therefore use of the above action limits us to even Nf . If, however, the fermions are integrated out, then the resulting effective action for the Aµ fields reads: Sef f = Nf Nf X 2 ln detM[A]. Aµ (x) − 2 8g xµ 2 (5) The effective action, though non-local, has an analytic dependence on Nf , and hence can be employed for odd or non-integer values.1 It can be simulated using a hybrid molecular dynamics algorithm, which evolves the {A} configuration through a fictitious time τ by a combination of microcanonical and Langevin dynamics. In the limit of timestep ∆τ → 0 the {A} are distributed according to the equilibrium ensemble. In this work we have implemented the R algorithm of Gottlieb et al [11]; in this case the systematic errors are O(N 2 ∆τ 2 ) [11, 12]. We found that on a 123 lattice a value ∆τ = 0.01 was sufficient, with a mean interval τ̄ = 1.0 between refreshments. Our results typically arise from averages over 500 units of τ . Fig. 1 shows a comparison between chiral condensate hχ̄χi data obtained for Nf = 4, for 1 Note that the form of (4) permits even-odd partitioning, which means that we can simulate integer powers of detM, rather than detM† M, with a local action. 3 0.16 0.14 0.12 σ 0.10 R algorithm HMC algorithm 0.08 0.06 0.01 0.02 0.03 0.04 0.05 m Figure 1: Chiral condensate (here denoted σ) vs. m on a 123 lattice for Nf = 4 with 1/g 2 = 1.0. various values of the bare mass m, using both the hybrid algorithm and a hybrid Monte Carlo algorithm [2], which in principle is free from systematic error. In Sec. 2.1 we will present results from simulations with variable Nf at a fixed coupling 1/g 2 = 1.0. The phase transition is located and found to fall at a value of Nf intermediate to the cases studied in earlier work [2, 3]. Fits to a power law equation of state ansatz yield critical exponents consistent with this picture. In Sec. 2.2 we present results obtained with variable g 2 at a fixed Nf = 5, the goal being to determine if the transition remains continuous, or whether there is evidence for coexisting phases signalling a first order transition, as observed for Nf = 6 [3]. Our conclusions, and the resulting phase diagram, are given in Sec. 3. 2 2.1 Numerical Simulations Fixed 1/g 2 = 1.0 Using the algorithm as described above we performed simulations with variable Nf on a 123 lattice at fixed 1/g 2 = 1.0, with bare mass m = 0.04, 0.03, 0.02. In Fig. 2 we plot the transverse and longitudinal susceptibilities χt,l , using the definitions given in 4 20.0 χt (directly computed) χt (computed with Ward identity) χl 15.0 χ 10.0 5.0 0.0 1.8 2.2 2.6 3.0 3.4 3.8 4.2 Nf Figure 2: Transverse and longitudinal susceptibilities vs. Nf on a 123 lattice with 1/g 2 = 1.0, m = 0.02. [2, 3]. As described there, χt can be calculated either from the integrated two-point pion propagator, or from hχ̄χi via the axial Ward identity: the two methods give compatible results, the latter being the less noisy. These quantities are related to the inverse square masses of respectively the pseudoscalar and scalar bound states. For Nf ≥ 3.6 they are roughly equal, indicating that the scalar and pion states are degenerate, and chiral symmetry is realised. For Nf ≤ 2.8, in contrast, χt ≫ χl , indicating that the pion has become light as expected for a pseudo-Goldstone mode in a phase of broken chiral symmetry. Fig. 2 thus offers evidence for a chiral symmetry breaking phase transition in the region Nf ≃ 3.0. Fig. 3 shows results for the order parameter hχ̄χi as a function of Nf , together with fits to a renormalisation-group inspired equation of state [2] of the form m = Athχ̄χiδ−1/β + Bhχ̄χiδ , (6) where δ, β are conventional critical exponents, and the parameter t expresses the distance from criticality. In previous work at fixed Nf , variable 1/g 2 we have defined t = 1/g 2 − 1/gc2 ; here we define t ≡ Nf − Nf ∗ , 5 (7) 0.4 m=0.02 m=0.03 m=0.04 m=0.00 0.3 σ 0.2 0.1 0.0 1.8 2.2 2.6 3.0 3.4 3.8 4.2 Nf Figure 3: Chiral condensate σ vs. Nf on a 123 lattice with 1/g 2 = 1.0, showing fits to the equation of state (6). The dashed line shows the fit in the chiral limit. where Nf ∗ denotes the value of Nf at the transition. Experience from previous studies [1, 2, 3] has shown that equations of state of the sort (6) provide adequate descriptions of the data in the vicinity of the transition, and that the resulting fitted exponents do not differ much from those extracted by more sophisticated finite volume scaling analyses. Also following previous work we fix β by the requirement δ − 1/β = 1, in order to stabilise the fitting procedure; we refer to the resulting four parameter fit as fit I. The results obtained by fitting data for Nf ∈ [3.0, 3.6] are given in Table 1. Also shown is a five parameter fit II over the same range for which the constraint on β was relaxed. A number of comments are in order. First, the constraint of fit I appears to be well satisfied by the data, since from fit II the fitted value of δ − 1/β = 0.97(8). Second, the fitted values of the exponents fall between those found for Nf = 2 and Nf = 4 from runs with variable 1/g 2 [1], being numerically closer to the Nf = 4 case: the results for these from comparable fits are summarised in Table 2. This supports the idea that gc2 (Nf ) is a critical line of fixed points along which critical exponents vary smoothly, with the numerical value of the exponents independent of the direction from within the (g 2, Nf ) plane from which the critical line is approached. Third, as a note of caution, it should be noted that the criterion for the range of Nf over which 6 Table 1: Results for 1/g 2 = 1.0 from fits to (6) on a 123 lattice. Parameter Fit I Fit II Nf ∗ 3.060(30) 3.060(44) A 0.264(7) 0.253(8) B 3.9(9) 4.4(1.3) δ 3.14(17) 3.22(24) β — 0.44(5) χ2 /d.o.f. 1.3 1.0 the fit was made was purely on the basis on minimising the χ2 per degree of freedom, as in [1, 2, 3]. Inspection of Fig. 3 suggests that the resulting fits may not describe the broken phase particularly well, and indeed fits which included smaller Nf values resulted in systematically larger values of δ. Our philosophy is to regard (6) merely as a phenomenological decription of the data; results from comparable fitting procedures applied to different models can thus be used to establish trends regardless of whether the true equation of state is ultimately of the form (6) or not. Table 2: Comparison with 123 lattice from ref. [1]. Parameter 1/gc2 (fit I) δ (fit I) β (fit II) 2.2 critical parameters for the Nf = 2 and 4 models from a Nf = 2 1.94(4) 2.68(16) 0.71(9) “Nf = 3.06” 1.00(1) 3.14(17) 0.44(5) Nf = 4 0.66(1) 3.43(19) 0.38(4) Fixed Nf = 5 With the confidence that simulations with a non-integer number of staggered flavors yield results consistent with our expectations for Nf ≃ 3, we then conducted a series of simulations on 123 lattices with similar parameters but this time with Nf = 5, a value which previous studies suggest is close to Nf c . Our results for hχ̄χi are shown in Fig. 4. We have included results from simulations with m = 0.01, although it is likely that these are badly finite-volume affected; however, they do reveal evidence of discontinuous behaviour for 1/g 2 ≃ 0.5. Fits to the equation of state (6), with the parameter t now given by 1/g 2 − 1/gc2, were attempted, but were less successful than those of the previous section, in general either resulting in an unacceptably large χ2 or not reaching convergence. In any case, the resulting fitted parameters are 7 0.25 m=0.02 m=0.03 m=0.04 m=0.01 0.20 0.15 σ 0.10 0.05 0.00 0.30 0.40 0.50 1/g 2 0.60 0.70 0.80 Figure 4: Chiral condensate σ as a function of 1/g 2 for Nf = 5 on a 123 lattice. strongly dependent on the fit range, and inclusion of the m = 0.01 points makes fit quality much worse. Our best results, coming from fitting data with m ∈ [0.04, 0.02], 1/g 2 ∈ [0.4, 0.6], are shown in Tab. 3. The fitted values of δ and β are consistent Table 3: Results for Nf = 5 from fits to (6) on a 123 lattice. Parameter Fit I Fit II 2 1/gc 0.437(4) 0.47(2) A 1.15(3) 0.5(2) B 190(40) 22(20) δ 5.6(1) 4.1(6) β — 0.29(5) 2 χ /d.o.f. 2.7 2.3 with the trends of Tab. 2, although it is interesting to note that the fitted value of δ−1/β = 0.6(2), indicating that in contrast to previous fits the data do not satisfy the usual constraint in this case. It could be argued that this disfavours the fit, since the combined assumptions of the degeneracy of scalar and pseudoscalar bound states in the chiral limit of the symmetric phase, plus the applicability of a power-law equation of state (6), naturally predict δ − 1/β = 1 [2, 13]. 8 0.08 0.06 σ 2 2 1/g =0.4 =0.45 =0.5 =0.52 =0.55 =0.6 =0.65 =0.7 0.04 0.02 0.00 0.00 0.05 0.10 0.15 0.20 0.25 0.30 m/σ Figure 5: Fisher plot of σ 2 vs. m/σ for Nf = 5 on a 123 lattice. To gain more insight we present the same data in the form of a Fisher plot (ie. hχ̄χi2 vs. m/hχ̄χi) in Fig. 5. This plot is devised to yield trajectories of constant 1/g 2 which intersect the vertical axis in the chiral limit in the broken phase, and the horizontal axis in the symmetric phase. Departures from the mean field indices δ = 3, β = 21 are revealed as curvature in the lines, and departures from δ − 1/β = 1 are revealed by variations in the sign of the curvature between the two phases. If we ignore the m = 0.01 points, which are those closest to the horizontal axis, then the plot suggests a critical 1/gc2 ≃ 0.5, with tentative evidence for δ − 1/β < 1. Inclusion of the low mass points, however, suggests an accumulation of the constant coupling trajectories around a line which if continued would intercept the horizontal axis. This is similar to the Fisher plot for 123 data for the model with Nf = 6, shown in fig. 6 of [2], and suggests similarities between the two cases. For Nf = 6, tentative evidence for tunnelling between coexisting vacua on simulations with m = 0.01 on a 163 lattice in the critical region was presented in [3], consistent with the chiral transition being first order. We tested this possibility by performing long Nf = 5 simulations on a 163 lattice for m = 0.01 and 1/g 2 = 0.45, 0.47 (these runs requiring ∆τ = 0.002): a time history for the latter is shown in Fig. 6. Whilst the fluctuations are large, with excursions of O(50%) about the mean, there is no convincing evidence for a two-state 9 0.25 0.20 σ 0.15 0.10 0.05 0 200 400 600 800 1000 simulation time Figure 6: Time history of σ for Nf = 5, m = 0.01 on a 163 lattice. signal. Double gaussian fits to the histograms of the binned data, which revealed coexisting states via twin peaks in [3], were this time statistically indistinguishable from simple single gaussians. In summary, whilst there seems to be clear evidence for a chiral transition with 1/g 2 = 0.50(5), the simulations performed to date are unable to determine the order. Fits based on the assumption of a continuous transition, whilst falling into the broad trend of existing data, are unsatisfactory and fail to satisfy the constraint δ −1/β = 1. A search for evidence of metastability consistent with a first order transition, on the other hand, proved negative. 3 The Phase Diagram In this Letter we have used a fresh simulation algorithm to explore non-integer numbers of staggered fermion flavors, and found results broadly consistent with earlier studies, namely that for Nf ≃ 3 the transition is continuous with critical exponents intermediate between those of Nf = 2 and Nf = 4 models, and that for Nf = 5 the order of the transition is difficult to determine, but shares features in common with the Nf = 6 transition, believed to be first order. This change in order of the 10 8 first order indeterminate order second order 6 Nf 4 2 0 0.0 0.5 1.0 1/g 2 1.5 2.0 2.5 Figure 7: Phase diagram for the three dimensional Thirring model; the chirally broken phase is to the lower left. The dashed line shows a fit to (9). transition with increasing Nf has also been observed in studies of QED4 [14]. The results also support the view that the partition function can be regarded as analytic in Nf , and that the non-locality of the fermionic action for Nf odd has no severe consequences. In this final section we collect together these new results with the best estimates obtained from previous finite volume scaling studies in [2, 3], to plot the phase diagram of the Thirring model on the (1/g 2, Nf ) plane in the chiral limit. The result is shown in Fig. 7. The shading of the points indicates the suspected order of the transition. The picture we have established is consistent with a critical line of fixed points, along which critical exponents vary smoothly, up to some Nf c ≃ 5, whereupon the transition becomes first order. It is work comparing our findings to the analytic approximations of Refs. [7] and [8]. Certainly the curvature of our critical line 1/gc2(Nf ) is consistent with both predictions (2) and (3). A detailed comparison with the former, which predicts a conformal fixed point at Nf = Nf c ≃ 4.32 is hampered by the difficulties in identifying the strong coupling limit in the lattice model, since there is an additive renormalisa2 2 tion relating 1/gCON T to 1/gLAT T due to non-conservation of the interaction current 11 in the lattice model [1, 2]. Therefore we are forced to identify Nf c with that value at which the transition becomes first order, at which point no continuum limit exists for the lattice model. Our estimate for Nf c is thus roughly consistent; however, the trends in the exponents extracted from the lattice studies, summarised in Table. 4 do not support a conformal fixed point. At such a point we expect the exponent δ to take the value 1, as for an asymptotically-free theory. The exponent η is related to the anomalous dimension of the composite operator χ̄χ by η = d − 2γχ̄χ [15], which in turn is predicted at the fixed point to be given by [8]: γχ̄χ = d−2 . 2 (8) At a conformal fixed point η should thus take the value 2. Finally, it is interesting to Nf 2 3.06(3) 4 5 Table 4: Critical exponents as a function of Nf . δ η 2.75(9) 0.60(2) 3.14(17) — 3.76(14) 0.26(4) 5.6(1)? — note that a fit of the form 1/gc2 ! π2 = A + B exp − Nf , 16 (9) motivated by Eq. (3) adjusted to allow for the renormalisation of 1/g 2 , can be made to the data shown in Fig. 7; the fit is particularly successful if the Nf = 2 point is excluded, in which case the resulting coefficients are A = 0.28(2), B = 4.7(2), with χ2 /d.o.f.=0.6. This suggests that the scenario of [7], in which chiral symmetry breaking is predicted for all Nf and hence that Nf c = ∞, cannot be excluded by the lattice studies to date. Acknowledgements This project was supported in part by the TMR-network “Finite temperature phase transitions in particle physics” EU-contract ERBFMRX-CT97-0122. We are greatly indebted to the CRT Computer Center of ENEL (PISA) for collaboration in the use of their CRAY YMP-2E. One of us (B.L.) would like to thank the Physics Department of Swansea for kindly hospitality. 12 References [1] L. Del Debbio and S.J. Hands, Phys. Lett. B373 (1996) 171. [2] L. Del Debbio, S.J. Hands and J.C. Mehegan, Nucl. Phys. B502 (1997) 269. [3] L. Del Debbio and S.J. Hands, hep-lat/9902014, to appear in Nucl. Phys. B. [4] M. Gomes, R.S. Mendes, R.F. Ribeiro and A.J. da Silva, Phys. Rev. D43 (1991) 3516. [5] S.J. Hands, Phys. Rev. D51 (1995) 5816. [6] T. Itoh, Y. Kim, M. Sugiura and K. Yamawaki, Prog. Theor. Phys. 93 (1995) 417. [7] D.K. Hong and S.H. Park, Phys. Rev. D 49 (1994) 5507. [8] M. Sugiura, Prog. Theor. Phys. 97 (1997) 311. [9] V.A. Miranskii and K. Yamawaki, Phys. Rev. D55 (1997) 5051. [10] C.J. Burden and A.N. Burkitt, Europhys. Lett. 3 (1987) 545. [11] S. Gottlieb, W.Liu, D. Toussaint, R.L. Renken and R.L. Sugar, Phys. Rev. D35 (1987) 2531. [12] S.J. Hands and J.B. Kogut, Nucl. Phys. B520 (1998) 382. [13] A. Kocić, J.B. Kogut and K.C. Wang, Nucl. Phys. B398 (1993) 405; M. Göckeler, R. Horsley, V. Linke, P.E.L. Rakow, G. Schierholz and H. Stüben, Nucl. Phys. B487 (1997) 313. [14] J.B. Kogut and K.C. Wang, Phys. Rev. D53 (1996) 1513. [15] S.J. Hands, A. Kocić and J.B. Kogut, Ann. Phys. 224 (1993) 29. 13