# 1 Introduction and Algorithm

## 1 Introduction and Algorithm

Recent numerical investigations  have revealed an interesting phase structure for the three dimensional Thirring model as a function of coupling $`g^2`$ and number of fermion species $`N_f`$. For sufficiently small $`N_f<N_{fc}`$, there is a continuous transition at $`g^2=g_c^2(N_f)`$ between a weak coupling phase in which chiral symmetry is realised in the limit $`m0`$, and a strong coupling one in which the symmetry is spontaneously broken. The critical indices characterising the transition are distinct for $`N_f=2`$ and $`N_f=4`$, suggesting that the continuum limits defined at the critical points define distinct interacting field theories. For $`N_f=6`$, however, simulations on a $`16^3`$ lattice with $`m=0.01`$ give tentative evidence for a first order chiral transition , implying both that no continuum limit exists in this case, and that $`4<N_{fc}<6`$.

Both the nature of the transition and the value of $`N_{fc}`$ are non-perturbative issues, inaccessible via either a standard perturbative expansion (which is non-renormalisable for $`d>2`$), or a $`1/N_f`$ expansion . There have, however, been analytic attempts to investigate the transition via the truncated Schwinger-Dyson (SD) equations . In the most systematic treatment , the SD equations are solved in ladder approximation in the the strong coupling limit; chiral symmetry breaking solutions are found for $`N_f<N_{fc}=128/3\pi ^24.32`$. Note that in a continuum approach the value of $`N_{fc}`$ manifests itself as the value of $`N_f`$ below which a non-trivial solution can be found, eg. by bifurcation theory. In a later paper , 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 $`N_f`$ read:

$`g^21:g_c^2`$ $`=`$ $`{\displaystyle \frac{2\pi ^2}{3}}N_f`$ (1)   
$`g^21:g_c^2`$ $`=`$ $`6\left[\mathrm{exp}\left({\displaystyle \frac{1}{\omega }}\left(2\pi 4\mathrm{tan}^1\omega \right)\right)1\right]`$   
with $`\omega ^2(N_f)={\displaystyle \frac{N_{fc}}{N_f}}1`$   
The essentially singular behaviour seen in (1) as $`N_fN_{fc}`$ is consistent with the existence of a conformal fixed point . Using a different sequence of truncations of the SD equations, however, Hong and Park  found chiral symmetry breaking solutions for all $`N_f`$, and predicted

$$g_c^2\mathrm{exp}\left(\frac{\pi ^2}{16}N_f\right),$$
(3)   
which is manifestly non-perturbative in $`1/N_f`$.

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 $`N_f`$ falling between the values $`N_f=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 $`N_{fc}`$. To start with, let us define the lattice action:

$`S`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{x\mu i}{}}\overline{\chi }_i(x)\eta _\mu (x)\left[(1+iA_\mu (x))\chi _i(x+\widehat{\mu })(1iA_\mu (x\widehat{\mu }))\chi _i(x\widehat{\mu })\right]`$   
$`+m{\displaystyle \underset{xi}{}}\overline{\chi }_i(x)\chi _i(x)+{\displaystyle \frac{N}{4g^2}}{\displaystyle \underset{x\mu }{}}A_\mu ^2(x)`$   
$``$ $`{\displaystyle \underset{xyi}{}}\overline{\chi }_i(x)_{ij}[A](x,y)\chi _j(y)+{\displaystyle \frac{N}{4g^2}}{\displaystyle \underset{x\mu }{}}A_\mu ^2(x),`$   
where the indices $`i,j`$ run over $`N`$ flavors of staggered lattice fermion. The $`A_\mu `$ are real auxiliary vector fields defined on the lattice links; gaussian integration over $`A_\mu `$ yields a form of the action with explicit four-fermion couplings  (other variants of the lattice-regularised Thirring model are discussed in ). Analysis of the spin-flavor content of staggered lattice fermions  reveals that in 3 dimensions the number of continuum four-component physical flavors is given by $`N_f=2N`$. Therefore use of the above action limits us to even $`N_f`$. If, however, the fermions are integrated out, then the resulting effective action for the $`A_\mu `$ fields reads:

$$S_{eff}=\frac{N_f}{8g^2}\underset{x\mu }{}A_\mu ^2(x)\frac{N_f}{2}\mathrm{ln}\text{det}[A].$$
(5)   
The effective action, though non-local, has an analytic dependence on $`N_f`$, and hence can be employed for odd or non-integer values.<sup>1</sup><sup>1</sup>1Note that the form of (1) permits even-odd partitioning, which means that we can simulate integer powers of $`\text{det}`$, rather than $`\text{det}^{}`$, with a local action. It can be simulated using a hybrid molecular dynamics algorithm, which evolves the $`\{A\}`$ configuration through a fictitious time $`\tau `$ by a combination of microcanonical and Langevin dynamics. In the limit of timestep $`\mathrm{\Delta }\tau 0`$ the $`\{A\}`$ are distributed according to the equilibrium ensemble. In this work we have implemented the R algorithm of Gottlieb et al ; in this case the systematic errors are $`O(N^2\mathrm{\Delta }\tau ^2)`$ . We found that on a $`12^3`$ lattice a value $`\mathrm{\Delta }\tau =0.01`$ was sufficient, with a mean interval $`\overline{\tau }=1.0`$ between refreshments. Our results typically arise from averages over 500 units of $`\tau `$. Fig. 1 shows a comparison between chiral condensate $`\overline{\chi }\chi `$ data obtained for $`N_f=4`$, for various values of the bare mass $`m`$, using both the hybrid algorithm and a hybrid Monte Carlo algorithm , which in principle is free from systematic error.

In Sec. 2.1 we will present results from simulations with variable $`N_f`$ at a fixed coupling $`1/g^2=1.0`$. The phase transition is located and found to fall at a value of $`N_f`$ intermediate to the cases studied in earlier work . 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 $`N_f=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 $`N_f=6`$ . Our conclusions, and the resulting phase diagram, are given in Sec. 3.

## 2 Numerical Simulations

### 2.1 Fixed $`1/g^2=1.0`$

Using the algorithm as described above we performed simulations with variable $`N_f`$ on a $`12^3`$ 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 $`\chi _{t,l}`$, using the definitions given in . As described there, $`\chi _t`$ can be calculated either from the integrated two-point pion propagator, or from $`\overline{\chi }\chi `$ 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 $`N_f3.6`$ they are roughly equal, indicating that the scalar and pion states are degenerate, and chiral symmetry is realised. For $`N_f2.8`$, in contrast, $`\chi _t\chi _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 $`N_f3.0`$.

Fig. 3 shows results for the order parameter $`\overline{\chi }\chi `$ as a function of $`N_f`$, together with fits to a renormalisation-group inspired equation of state  of the form

$$m=At\overline{\chi }\chi ^{\delta 1/\beta }+B\overline{\chi }\chi ^\delta ,$$
(6)   
where $`\delta ,\beta `$ are conventional critical exponents, and the parameter $`t`$ expresses the distance from criticality. In previous work at fixed $`N_f`$, variable $`1/g^2`$ we have defined $`t=1/g^21/g_c^2`$; here we define

$$tN_fN_f,$$
(7)   
where $`N_f`$ denotes the value of $`N_f`$ at the transition. Experience from previous studies  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 $`\beta `$ by the requirement $`\delta 1/\beta =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 $`N_f[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 $`\beta `$ 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 $`\delta 1/\beta =0.97(8)`$. Second, the fitted values of the exponents fall between those found for $`N_f=2`$ and $`N_f=4`$ from runs with variable $`1/g^2`$ , being numerically closer to the $`N_f=4`$ case: the results for these from comparable fits are summarised in Table 2. This supports the idea that $`g_c^2(N_f)`$ 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,N_f)`$ 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 $`N_f`$ over which the fit was made was purely on the basis on minimising the $`\chi ^2`$ per degree of freedom, as in . Inspection of Fig. 3 suggests that the resulting fits may not describe the broken phase particularly well, and indeed fits which included smaller $`N_f`$ values resulted in systematically larger values of $`\delta `$. 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.

### 2.2 Fixed $`N_f=5`$

With the confidence that simulations with a non-integer number of staggered flavors yield results consistent with our expectations for $`N_f3`$, we then conducted a series of simulations on $`12^3`$ lattices with similar parameters but this time with $`N_f=5`$, a value which previous studies suggest is close to $`N_{fc}`$. Our results for $`\overline{\chi }\chi `$ 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^20.5`$. Fits to the equation of state (6), with the parameter $`t`$ now given by $`1/g^21/g_c^2`$, were attempted, but were less successful than those of the previous section, in general either resulting in an unacceptably large $`\chi ^2`$ or not reaching convergence. In any case, the resulting fitted parameters are 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 $`\delta `$ and $`\beta `$ are consistent with the trends of Tab. 2, although it is interesting to note that the fitted value of $`\delta 1/\beta =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 $`\delta 1/\beta =1`$ .

To gain more insight we present the same data in the form of a Fisher plot (ie. $`\overline{\chi }\chi ^2`$ vs. $`m/\overline{\chi }\chi `$) 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 $`\delta =3,\beta =\frac{1}{2}`$ are revealed as curvature in the lines, and departures from $`\delta 1/\beta =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/g_c^20.5`$, with tentative evidence for $`\delta 1/\beta <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 $`12^3`$ data for the model with $`N_f=6`$, shown in fig. 6 of , and suggests similarities between the two cases. For $`N_f=6`$, tentative evidence for tunnelling between coexisting vacua on simulations with $`m=0.01`$ on a $`16^3`$ lattice in the critical region was presented in , consistent with the chiral transition being first order. We tested this possibility by performing long $`N_f=5`$ simulations on a $`16^3`$ lattice for $`m=0.01`$ and $`1/g^2=0.45`$, 0.47 (these runs requiring $`\mathrm{\Delta }\tau =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 signal. Double gaussian fits to the histograms of the binned data, which revealed coexisting states via twin peaks in , 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 $`\delta 1/\beta =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 $`N_f3`$ the transition is continuous with critical exponents intermediate between those of $`N_f=2`$ and $`N_f=4`$ models, and that for $`N_f=5`$ the order of the transition is difficult to determine, but shares features in common with the $`N_f=6`$ transition, believed to be first order. This change in order of the transition with increasing $`N_f`$ has also been observed in studies of $`\text{QED}_4`$ . The results also support the view that the partition function can be regarded as analytic in $`N_f`$, and that the non-locality of the fermionic action for $`N_f`$ 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 , to plot the phase diagram of the Thirring model on the $`(1/g^2,N_f)`$ 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 $`N_{fc}5`$, whereupon the transition becomes first order.

It is work comparing our findings to the analytic approximations of Refs.  and . Certainly the curvature of our critical line $`1/g_c^2(N_f)`$ is consistent with both predictions (1) and (3). A detailed comparison with the former, which predicts a conformal fixed point at $`N_f=N_{fc}4.32`$ is hampered by the difficulties in identifying the strong coupling limit in the lattice model, since there is an additive renormalisation relating $`1/g_{CONT}^2`$ to $`1/g_{LATT}^2`$ due to non-conservation of the interaction current in the lattice model . Therefore we are forced to identify $`N_{fc}`$ with that value at which the transition becomes first order, at which point no continuum limit exists for the lattice model. Our estimate for $`N_{fc}`$ 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 $`\delta `$ to take the value 1, as for an asymptotically-free theory. The exponent $`\eta `$ is related to the anomalous dimension of the composite operator $`\overline{\chi }\chi `$ by $`\eta =d2\gamma _{\overline{\chi }\chi }`$ , which in turn is predicted at the fixed point to be given by :

$$\gamma _{\overline{\chi }\chi }=\frac{d2}{2}.$$
(8)   
At a conformal fixed point $`\eta `$ should thus take the value 2.

Finally, it is interesting to note that a fit of the form

$$1/g_c^2=A+B\mathrm{exp}\left(\frac{\pi ^2}{16}N_f\right),$$
(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 $`N_f=2`$ point is excluded, in which case the resulting coefficients are $`A=0.28(2)`$, $`B=4.7(2)`$, with $`\chi ^2`$/d.o.f.=0.6. This suggests that the scenario of , in which chiral symmetry breaking is predicted for all $`N_f`$ and hence that $`N_{fc}=\mathrm{}`$, 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.