{'fraction_non_alphanumeric': 0.03517830097755748, 'fraction_numerical': 0.01540341456698334, 'mean_word_length': 4.297684172137126, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Three Dimensional Thirring Model for Small N_f

We formulate the three dimensional Thirring model on a spacetime lattice and study it for various even numbers of fermion flavors N_f by Monte Carlo simulation. We find clear evidence for spontaneous chiral symmetry breaking at strong coupling, contradicting the predictions of the 1/N_f expansion. The critical point appears to correspond to an ultra-violet fixed point of the renormalisation group; a fit to a RG-inspired equation of state in the vicinity of the fixed point yields distinct critical exponents for N_f=2 and N_f=4, while no fit is found for N_f=6, suggesting there is a critical number N_fc<6 beyond which no chiral symmetry breaking occurs. The spectrum of the N_f=2 theory is studied; the states examined vary sharply but continuously across the transition.


Introduction
The Thirring model is a quantum field theory of fermions interacting via a currentcurrent contact term, described in three dimensional continuum Euclidean spacetime by the following Lagrangian: (1.1) Here ψ andψ are four-component spinors, m is a bare mass, and i runs over N f fermion species. After the introduction of an auxiliary vector field, the above Lagrangian can be rewritten as: which now contains the fermionic bilinear term of the usual QED lagrangian but with a mass term for the vector field, which spoils the gauge-invariance of L aux . This is no surprise, since the original lagrangian did not show any local symmetry. The global symmetries are a U(N f ) one, corresponding to rotations in flavor space, and, for m → 0, a U(N f ) A chiral symmetry. Since the coupling constant g has mass dimension [g] = (1− d 2 ) = − 1 2 , the usual perturbative expansion is non-renormalisable. However perturbation theory is not the end of the story and, as suggested several years ago [1,2], different approaches can yield a well-defined continuum limit, corresponding to a UV-stable fixed point of the renormalisation group (RG).
The Gross-Neveu model (GN), which has been studied both numerically and in the 1/N f expansion [3,4,5], exhibits this behaviour. The model turns out to be renormalisable to all orders in 1/N f , with a UV-stable RG fixed point at g U V = 0 and chiral symmetry breaking for g > g U V ; only for g ≃ g U V can the ratio of the physical scale (either a particle mass or a scattering cross-section) to the UV cutoff be made small. The theory defined at this point describes the IR dynamics of a linear sigma model with a Yukawa coupling to fermions [6], which is super-renormalisable about the origin: For a critical value of v in the IR regime the kinetic and quartic terms of the σ-field become irrelevant operators and thus the system equivalent to the GN model. The critical behaviour for the GN model can be summarised by the behaviour of the βfunction, as shown in Fig. 1. In two dimensions, the theory is asymptotically free, and the β-function exhibits the usual UV fixed point at the origin. For 2 < d < 4, the fixed point is shifted away from the origin, separating a strong coupling regime, where dynamical chiral symmetry breaking is possible, from a weak coupling phase where the theory is chirally symmetric [7]. The situation for the Thirring model is less clear. At both leading and next-toleading order in 1/N f the renormalisability of the massless model has been established [8]. For d < 4 the vacuum polarisation is actually UV-finite so long as the regularisation respects current conservation. As a consequence a continuum limit may be defined for any value of the coupling g, which would control, e.g., the ratio of physical fermion mass to vector bound state mass. In other words, the β-function is vanishing for any value of g and the interaction (ψ i γ µ ψ i ) 2 is a marginal operator, whereas the interaction for the GN model is relevant [9]. Another prediction of the 1/N f expansion is the vanishing of the vector bound state mass in the strong coupling limit mg → ∞. A massless particle mediating an interaction between conserved currents suggests an equivalence between the UV behaviour of Thirring model and the IR one of QED 3 , which would thus be the analogue of the equivalence between the GN and Yukawa models. Such a conjecture is supported by the fact that the O(1/N f ) corrections in the two models appear to coincide [8,10]. Moreover in the 1/N f expansion, as in perturbative gauge theory, no diagram ever appears resulting in the dynamical breaking of chiral symmetry leading to spontaneous mass generation.
The prediction of no spontaneous mass generation is contradicted, however, by results from a completely non-perturbative approach, namely solution of Schwinger-Dyson (SD) equations, which predicts a non-vanishing chiral condensate ψ ψ in the limit m → 0 [2,11,12]. In this approach a non-trivial solution for the dressed fermion propagator S(p) = (A(p)ip / + Σ(p)) −1 is sought, i.e. one in which the self-energy Σ(p) and hence ψ ψ are non-vanishing in the chiral limit m → 0. Unfortunately, the SD equations can only be solved by truncating them in a somewhat arbitrary fashion. The usual approximation [2,12], is to assume that the vector propagator is given by its leading-order form for m = 0 in the 1/N f expansion, viz.
and that the fermion-vector vertex function is well-approximated by the bare vertex (the so-called "planar" or "ladder" approximation): The most systematic treatment has been given by Itoh et al [12], who note the equivalence of (1.2) to a gauge-fixed form of a gauged fermion-scalar model. This allows the specification of an alternative non-local gauge-fixing condition, which in turn entails A(p) ≡ 1, with no effect, of course, on gauge-invariant quantities such as the physical mass or ψ ψ . The resulting SD equations can be solved exactly in the limit g 2 → ∞; a non-trivial solution for Σ(p) exists for This is identical to the critical N f c predicted for non-trivial IR behaviour in QED 3 [13]. However, since the integral equations require the introduction of a UV cutoff Λ, a feature of this solution is that the induced physical scale µ depends on N f in an essentially singular way: This implies that in the strong coupling limit a continuum limit only exists as N f → N f c . The behaviour (1.7) is an example of an infinite order, or "conformal" phase transition, originally discussed in a particle physics context by Miranskii and collaborators in quenched QED 4 [14], and recently discussed in more generality [15]. One prediction of [15] is that although the order parameter ψ ψ is continuous at the transition, there is a discontinuity in the spectrum of light excitations, with no light scalar resonances on the symmetric side of the transition. In the broken phase, of course, Goldstone's theorem predicts light pions. This behaviour has recently been exhibited in an investigation of four dimensional SU(N c ) gauge theory with an intermediate number N f of fermion flavors [16]: a solution of the Bethe-Salpeter equation in ladder approximation found no low-momentum pole in the symmetric phase.
Unfortunately no analytic solution to the SD equation for the Thirring model exists for g 2 < ∞; however using different techniques Kondo [17] has argued that there is a critical line N f c (g 2 ) in the (g 2 , N f ) plane, which is a smooth invertible function. Therefore for integer N f < N f c one might conjecture the form [18] corresponding to a symmetry restoring transition at some critical point g 2 = g 2 c (N f ). Presumably in this scenario the coupling g has become relevant: there may exist a novel strongly-coupled continuum limit at the critical point not described by the 1/N f expansion, which is simultaneously the UV limit of the Thirring model and the IR limit of QED 3 . Since there is no small dimensionless parameter in play, this fixed point would be inherently non-perturbative. It is worth mentioning in passing that a novel IR fixed point in QED 3 has been proposed by Aitchison and Mavromatos in connection with non-Fermi liquid behaviour in the normal phase of superconducting cuprates [19].
Although this scenario of a non-trivial fixed point for N f < N f c is attractive, there are good reasons to be cautious. Using a different sequence of truncations Hong and Park have found chiral symmetry breaking for all N f [11], with a result which is also non-analytic in 1/N f and hence beyond the reach of the 1/N f perturbative series. Moreover in QED 3 , studies beyond the planar approximation, using improved ansätze for the vertex Γ µ , suggest that the condition A(p) ≡ 1 is unphysical, and that chiral symmetry is spontaneously broken for all N f [20]. In the current context this would imply g 2 c < ∞ for all N f . Because of the systematic uncertainty attached to the SD approach, and because of the intrinsic interest in exploring the non-perturbative behaviour of a fermionic model, we have studied the three dimensional Thirring model by Monte Carlo simulation on a spacetime lattice. Our initial results appeared in [18], where we claimed that critical points with distinct critical indices exist for the models with N f = 2, 4, but that evidence for a critical theory with N f = 6 was not found. In this paper we summarise those results, and present further data which enable us to refine our description of the critical behaviour of the N f = 2 model. In Sec. 2 we give a fairly lengthy exposition of the lattice formulation, showing that the lattice method is itself not without systematic uncertainties. Although we use staggered lattice fermions which have a continuous remnant of the U(N f ) A axial symmetry, a transcription to a basis in which the Fermi fields have explicit spin and flavor indices reveals that there are many other interactions present at tree level besides the desired current-current form. Perhaps a robust counter-argument to this problem is that for a strongly-coupled fixed point a non-perturbative regularisation is always necessary, and therefore "what you see is what you get", i.e. the fixed point theory defines itself via its critical exponents and spectral quantities. In addition, the lattice form of the current is not exactly conserved, resulting in an additive renormalisation of the inverse coupling 1/g 2 . This makes it difficult in practice to identify the strong coupling limit g 2 → ∞. A description of the simulation, and a list of the measured observables, namely the condensate, susceptibilities and propagators, follows.
In Sec. 3 we present results for the equation of state, i.e. m = f ( ψ ψ , g). Since it is impossible to work directly in the chiral limit m = 0, we find it convenient to fit to a phenomenological form inspired by the standard power-law scaling form used in the theory of ferromagnetism, with its attendant critical exponents β, γ, δ and ν, together with the additional constraints of hyperscaling, and, from the ladder approach to the gauged Nambu -Jona-Lasinio model [21], that δ − 1/β = γ ≡ 1. We justify this approach not so much by any theoretical argument but because its shortcomings if any should manifest themselves as variations in the fitted values of the exponents as data is taken from systems on larger volumes and closer to the chiral limit. In [18] the fitted exponents for the N f = 2, 4 models were distinct, a feature not seen in similar studies of QED 4 [22,23]. Here we supplement the data of [18] for N f = 2 with results from a larger lattice (16 3 ) and a smaller bare mass m, and refine the fit by combining data from several lattice volumes in a finite size scaling analysis. We find critical exponents and g 2 c completely consistent with those of our earlier studies, which strengthens the case for a critical point described by power-law scaling. We also have some new results for the N f = 6 model on 12 3 systems, supporting our previous claim that there is no transition in this case, and hence that N f c < 6 for the Thirring model.
In Sec. 4 we present exploratory results for susceptibilities, and for fermion and bound-state masses. Our lattice size of 16 3 is unfortunately not long enough in the temporal direction to permit real accuracy; nonetheless qualitative features of the level ordering in either phase do emerge. The susceptibilities reflect the phase transition very nicely, with longitudinal and transverse susceptibilities χ l and χ t varying from near-degeneracy in the symmetric phase (thus indirectly supporting the hypothesis δ − 1/β = 1) to χ t ≫ χ l in the broken phase. The effects of dynamical fermion loops in the quantum vacuum can also be estimated, and appear to become dominant in the broken phase. The fermion, pseudoscalar and scalar masses all behave as expected, the pion mass staying small and roughly constant, while fermion and scalar masses change sharply but continuously across the transition. There is no sign of the discontinuities described in [15]. Finally studies in the vector-like channel reveal light states both with vector and axial vector quantum numbers in the symmetric phase. This is interesting because no light axial vector is predicted in the continuum 1/N f expansion at leading order, as we show in Sec. 4. Both states become difficult to observe in the broken phase. Our conclusions are presented in Sec. 5.

Lattice formulation
We want to discretise the action of Eq. (1.1) on a hypercubic Euclidean lattice. In order to recover a partial U(N) A axial symmetry in the massless limit, we use N flavors of staggered fermions. The lattice action then reads: where χ,χ are the staggered fermion fields, η µ the Kawamoto-Smit phases, m is the bare fermion mass, the flavor indices i, j run from 1 to N and the vector character of the interaction is encoded in the geometrical structure of the four fermion term. To express the action as a bilinear in χ,χ, which makes it amenable to simulation, we introduce a bosonic auxiliary field as in (1.2). In this paper we choose a form which we will refer to as non-compact: where we have introduced M[A, m] to denote the fermionic bilinear, which depends on both the auxiliary field and the bare mass. It is straightforward to integrate over the auxiliary field and recover Eq. (2.1).
In the chiral limit m → 0 both the original action (2.1) and the bosonised form (2.2) are invariant under a global axial rotation, which for N = 1 is written where the phase ε(x) = (−1) x 1 +x 2 +x 3 . The generalisation to U(N) is straightforward.
We will refer to this as the chiral symmetry of the model. In the m → 0 limit for certain values of the coupling g 2 > g 2 c (N) it is spontaneously broken, signalled by the appearance of a condensate χχ = 0.
Some comments on the lattice formulation are appropriate. Eq. (2.2) in the long wavelength limit (i.e. assuming a∂ψ ≪ ψ) is clearly closely related to the continuum auxiliary form (1.2), and as such is suitable for a perturbative expansion in small-mode fluctuations of the A µ fields, corresponding to the continuum 1/N f expansion. This expansion is well-behaved because to leading order the vector propagator D µν (x) = A µ (0)A ν (x) receives a contribution not just from the quadratic term A 2 µ , but also from the vacuum polarisation loop Π µν , which is finite in d < 4. We find in the deep Euclidean region [8] where P µν is the transverse projector Because the interaction currentψγ µ ψ is conserved, the longitudinal component of D µν has no physical significance [1][12] [17]. The 1/k behaviour of D µν in effect suppresses short wavelength excitations of the auxiliary field, and justifies the long-wavelength expansion a posteriori. It means we might hope the form (2.2) has the correct continuum limit. This assumption may not be justified, however, due to the chiral phase transition at finite g 2 . Since the 1/N f expansion predicts no phase transition, we know that there will be at least some part of coupling constant space where the above arguments cannot hold. In this case it is not quite so clear that (2.2) coincides with the continuum form (1.1). To see why there might be a problem, we must consider the continuum flavor interpretation of the lattice action (2.1). First define a unitary transformation to fields q,q [24]: where y denotes a site on a lattice of spacing 2a and P is a 3-vector with entries either 0 or 1 ranging over corners of the elementary cube associated with y, so that each side of the original lattice corresponds to a unique choice of (P ; y). The 2 × 2 matrices Γ P and B P are defined by where τ µ are the Pauli matrices. It is then possible to recast the bilinear term of (2.1) as follows; The direct product is between 4 × 4 Dirac matrices, acting on spinor degrees of freedom, defined by q → exp(iα(γ 5 ⊗ 1 1))q ;q →q exp(iα(γ 5 ⊗ 1 1)). (2.11) The condensate χχ in the new basis reads q(1 1 ⊗ 1 1)q . The first and third terms of (2.8) are clearly of the same form as the continuum action for free massive fermions, whereas the second term, which is flavor non-singlet and Lorentz non-covariant, is formally O(a) and hence hopefully irrelevant. The problem comes when we consider the interaction term. In the auxiliary field formalism, only in the long-wavelength limit would we expect the interaction to assume its expected continuum form A µqi (γ µ ⊗ 1 1 2 )q i [25]. When A is integrated out to recover the original form (2.1) though, all momentum modes are summed over with equal weight. In fact, in four-fermi form the interaction contains many extra terms 2 : where the difference operators f and δ are defined by The terms containing δ are formally O(a) and hence probably irrelevant, but no such argument can be used for the remainder, which, in addition to the desired form (q i (γ µ ⊗ 1 1)q i ) 2 , contains extra flavor non-singlet and non-covariant terms. If the 1/N f expansion were valid, then one could argue that only interactions between conserved currents would generate long-range correlations between the A fields, and that all other interactions would become irrelevant in the continuum limit -a similar effect appears to be the case in other two [27,28] and three [3] dimensional four-fermi models where a renormalisable expansion is available. Even in this case there might also be contributions in the chiral limit from the interaction (q i (γ µ γ 5 ⊗ 1 1)q i ) 2 . However, as argued above, in the present case we know the 1/N f expansion fails qualitatively in the chirally broken phase. It therefore becomes a matter of experiment, by studying the current-current and AA correlations in various channels numerically, to determine the continuum limit of the action (2.1). A second issue we must address is that despite the above comments there is an important difference between the continuum 1/N f expansion and the 1/N expansion applied to the non-compact lattice action (2.2). This concerns vector current conservation. In the continuum model a regularisation such as Pauli-Villars, which conserves the interaction current, must be assumed in order to ensure that the superficial O(Λ) contribution to vacuum polarisation vanishes, leaving a finite result [1,2,8]. Now, in non-compact lattice QED 3 , the vacuum polarisation is again UV finite, but this time because of a cancellation between O(a −1 ) contributions from two separate diagrams ( Fig. 2 -see eg. [29]). The second diagram appears because the gauge-invariant interaction term in QED isχ exp(ieA µ )χ, and hence contains n-photon-electron-positron vertices. In the action (2.2), on the other hand, the interaction isχiA µ χ, so only the left-hand diagram is present. Therefore cancellation of O(a −1 ) divergences does not occur. We find for the vacuum polarisation tensor to leading order in 1/N: and, recalling a factor of two from (2.10) This means that for the inverse auxiliary propagator (setting a = 1) If we compare this with the continuum result we see that the O(a −1 ) divergence can be absorbed by a wavefunction renormalisation together with a coupling constant renormalisation: . (2.20) Now, the physics described by continuum 1/N f perturbation theory occurs for the range of couplings g 2 R ∈ [0, ∞), i.e. for g 2 ∈ [0, g 2 lim ), where to leading order We therefore expect discontinuous behaviour at bare coupling g 2 ≃ g 2 lim , since at this point the auxiliary propagator is not defined (i.e. D −1 µν in (2.17) has no inverse). For g 2 > g 2 lim , since Π CON T → 0 in the limit k 2 → 0, D µν becomes negative, suggesting that the model is not unitary in this region.
In Fig. 3 we plot σ = χχ vs. 1/g 2 for N = 1, 2, 3 and m = 0.10 on a 8 3 lattice [18]. The three models have apparently coincident condensates in the strong coupling region 1/g 2 ≤ 0.3, which appear to be tending to zero in the g 2 → ∞ limit. It is tempting to associate this behaviour with g 2 > g 2 lim , but the correspondence with the value (2.21) is not good. It may well be that the numerical value of the right-hand diagram of Fig. 2 is considerably reduced in a chirally broken vacuum; in any case we do not expect the 1/N expansion to be accurate here for the reasons described above.
The ill-defined behaviour of the g 2 → ∞ limit of the non-compact action (2.2) presents no problem for most of the simulation results presented in this paper, which focus on the critical region of the N = 1 model around 1/g 2 c ≃ 2.0, and of the N = 2 model with 1/g 2 c ≃ 0.6. However, we expect g 2 c to increase as N increases [11,12,17,18], so that as models with larger N are explored, eventually the relation between g 2 and g 2 R will become important, particularly if lattice studies are to resolve the issue of the existence and value of N f c . It is therefore worth making remarks on two other possible forms for the lattice Thirring model.
For N = 1 there is an alternative but equivalent way of introducing the auxiliary field, which we will call the compact form: where the link variables θ µ (x) are freely integrated over [30,31]. It is straightforward to perform this integration to yield a term ∝ I 0 ( −2g 2χ (x)χ(x +μ)χ(x +μ)χ(x)), with I 0 a modified Bessel function. If we define I 0 by its power series expansion, then the Grassmann nature of the χ,χ truncates the expansion and forces equivalence with the partition function based on (2.1). We have performed some trial simulations using (2.22), and found identical results to those of the non-compact form (2.2) within statistical accuracy. It is interesting to note in passing that using a similar method one can show the equivalence of (2.1) with the strong gauge-coupling (i.e. β = 0) limit of the gauge-fermion-scalar (χUφ) model recently studied in connection with non-trivial continuum limits in 3 and 4 dimensions [32]. The bilinear action (2.22) can be extended to N > 1 simply by introducing explicit flavor indices i on χ,χ. Now, however, terms of the form (χ(x)χ(x +μ)χ(x +μ)χ(x)) n with n > 1 are not suppressed by the Grassmann nature of χ: for N flavors of staggered fermion there are thus n-point couplings in the action with n ≤ 4N. Since ultimately we wish to understand and contrast the behaviour of the Thirring model for various values of N, we prefer the non-compact form (2.2), where the interaction is of the same form for all N. It may well prove to be the case that in the RG sense the higher-point couplings are irrelevant, and hence (2.2) and (2.22) have identical continuum limits. Both variants deserve study.
Finally we mention a recent numerical study of a gauge-invariant version of the three dimensional Thirring model by Kim and Kim [33]; the action is with φ a gauge-covariant scalar field. This action has a local gauge invariance by construction, corresponding to the hidden local symmetry noted in Refs. [12,17]. In unitary gauge φ = 1 this action resembles (2.2) in the long-wavelength limit, but now has a periodic dependence on the auxiliary field. The main result of [33] is that the N f = 6 model appears qualitatively different from N f = 2, once again supporting a value of N f c < 6.
It is an open question whether the action (2.23) lies in the same universality class as (2.1). We simulated dynamical fermions using a hybrid Monte Carlo (HMC) algorithm. In all cases we used a random trajectory length drawn from a Poisson distribution with mean 0.9, and adjusted the timestep to maintain an acceptance rate of 70% or greater. Typically we performed measurements on configurations separated by between 2 and 5 trajectories. Our estimated autocorrelation times measured on sequences of configurations varied from ∼ 1 to ∼ 3 near the critical coupling; all our quoted errors take autocorrelations into account. The convergence criterion we used in the conjugate gradient matrix inversion was that the modulus of the residual vec-tor was 10 −6 per lattice site during guidance and 10 −9 per site for Metropolis and measurement. Unlike the case of the GN model [3], we were unable to boost acceptance rates by tuning the parameters of the guidance Hamiltonian. However, the vectorlike nature of the interaction allows the simulation of the model for any value of the number of lattice species N, corresponding to a number of physical fermions N f = 2N, via an even-odd partitioning. For the very same reason, the only diagonal term in the fermion matrix comes from the mass term in the action. As m gets small, we find that the fermion matrix becomes very difficult to invert. Simulations with small m are therefore extremely time-consuming and represent a major problem for the extrapolation of the data to the chiral limit. Most of the results reported in this paper use between 200 and 500 measurements at each value of the coupling and bare mass. We have used periodic boundary conditions for the auxiliary field, and for the fermions used periodic in the spacelike directions and antiperiodic in the timelike direction. Most of the new results in the paper are obtained on 16 3 lattices; however we also did some runs on an asymmetric 8 2 ×16 lattice, to give some estimate of finite volume effects on correlation function measurements.
In our simulations we measured the following observables: (i) The chiral condensate where M is the fermion kinetic operator introduced in eq.(2.2). In Sec. 3 we will discuss numerical fits to the model's equation of state using these measurements. To enable alternative fitting procedures we list all our measurements made on symmetrical (i.e. L 3 ) lattices for N f = 2 in Tabs. 1 and 2, for N f = 4 in Tab. 3 and for N f = 6 in Tab. 4.
(ii) The longitudinal susceptibility where we have denoted the components with disconnected and connected fermion lines respectively as "singlet" χ ls and "non-singlet" χ lns in analogy with QCD flavor dynamics.
(iii) The transverse susceptibility where we have used the fact that χεχ vanishes, and the last line follows from a Ward identity.
(iv) The "vector susceptibility" where V µ is the interaction current: Due to the inhomogeneous boundary conditions, we only averaged this quantity over the spacelike directions. All susceptibilities, of course, are identical to the bound state propagator in the appropriate channel evaluated at zero momentum. The non-singlet components were therefore easily evaluated from the bound state correlators discussed below. The singlet parts, involving traces over a single power of M −1 , were evaluated using gaussian noise vectors. Since the susceptibility signals are intrinsically noisy, we used 5 noise vectors per configuration, thus obtaining 5 independent estimates of χχ , and 5 × 4/2 = 10 independent susceptibility estimates.
In addition to these operators we have also performed spectroscopy by evaluating timeslice propagators. We examined the following channels: where note that the signal is improved if only sites an even number of translations in the transverse directions are included in the sum.
(vii) One-link composite operators i.e. the two-point correlator between interaction currents V µ (x), corresponding to the non-singlet component of the vector susceptibility (2.27).
Of course, before assigning names to the various channels we should examine the spin/flavor structure of the fermion bilinears using the transformation (2.6). Since the staggered fermion action is invariant only under translations by an even number of lattice spacings, in general, as t → ∞, two-point correlators behave as i.e. displaying signals in both "direct" and "alternating" channels, characterised by distinct masses M d and M a respectively. The spin/flavor assignments we find are shown in Tab. 5. It is interesting to note that in three Euclidean dimensions the local rho meson operator, which we have referred to as "local vector", does not in fact have vector quantum numbers.
In the calculation of timeslice propagators we evaluated the inverse matrix M −1 using a single local source on each configuration, chosen at a site randomly displaced from the origin by an even number of lattice spacings in each direction.

The equation of state from RG equations
Spontaneous breaking of chiral symmetry is signalled by a non-vanishing chiral condensate as m → 0. As already pointed out in the previous section, it is extremely difficult to run the HMC for very small values of the bare mass. Therefore we are bound to work with a restricted set of values for m, making an extrapolation to m = 0 at a fixed value of g rather unreliable.
In order to find numerical evidence of the RG structure described in the introduction, the scaling properties of the theory at the UV fixed point are going to be used to find an equation of state allowing a global fit of the data. The equation of state relates the order parameter for the broken symmetry (in our case the chiral condensate) to the value of the external symmetry breaking field (the bare mass m) as a function of the external parameters (coupling constant, finite size of the lattice), thus providing a powerful tool for investigating the behaviour of the order parameter near the phase transition in the (g, m)-plane.

Scaling properties at fixed lattice size
Let us briefly recall how an equation of state can be obtained from the RG equations (RGE) (see eg. [34]). Consider a theory of a bosonic field φ regularised with some cut-off Λ, coupled to an external field H and let M = φ be the order parameter. Then the effective potential can be expanded in powers of M: where Γ (n) are the 1PI Green functions at zero external momenta, t is the reduced coupling (see (3.10)) and Λ the cut-off. The external field, for given magnetisation, can be written as: The above equation allows the RGE for H to be written: where: and Z is the wave-function renormalisation. The solution of the RGE is obtained by using the method of characteristics. In a neighbourhood of the fixed point, the solution can be combined with dimensional analysis in order to rescale the cut-off and obtain the general result: where F is a universal scaling function. By setting H = 0 in eq. (3.6), the critical behaviour of the order parameter when the external field is switched off is recovered: while, for t = 0, F (0) is a constant and hence: showing clearly that β and δ are the usual critical exponents introduced in the context of phase transitions. If the critical exponents are related to the existence of a UV fixed point, as we are assuming in this section, they must obey the hyper-scaling relations, obtained under the assumption that there is a single physical length scale near the continuum limit: where η is the anomalous dimension and ν is the critical exponent which characterises the divergence of the correlation length as t → 0. For the Thirring model the reduced coupling t is identified with: and the symmetry breaking field with the bare mass m, while the order parameter is the chiral condensate χχ defined in the previous section. A Taylor expansion for small t reduces Eq. (3.6) to our first equation of state: If the critical behaviour is described by mean-field theory, the critical exponents in Eq. (3.6) are expected to be δ = 3 and β = 1/2, leading to: In mean-field theory the square of the chiral condensate is a linear function of the ratio m/ χχ . A plot of χχ 2 vs. m/ χχ is called a Fisher plot. In such a plot a positive value of the intercept corresponds to a non-vanishing condensate for m = 0, i.e. to chiral symmetry breaking, while the value of the intercept will be exactly zero at the critical coupling. Eq. (3.11) permits a five parameter fit of the numerical data for the chiral condensate. In order to reduce the number of free parameters in the fit, one can set δ − 1/β = 1 and keep δ as the only exponent to be determined from the data. This relation originally arose from solution of the SD equations for the gauged Nambu -Jona-Lasinio model in four dimensions in ladder approximation [21]; here we simply use it as a plausible hypothesis to be tested against the data. In what follows we refer to the 5 and 4 parameter fits as fits I and II respectively. The outcome of fit I shows that the above relation between δ and β is satisfied wihin errors.
Since both formulae rely on a Taylor expansion around t = 0, they are only expected to describe the behaviour of the chiral condensate in a neighbourhood of the critical coupling, so that only a reduced set of values of g can be fitted by our equation of state. The number of points included in the fit has been chosen in such a way to minimise the χ 2 value.
In Figs. 4 and 5 we show the Fisher plots for N f = 2, 4 on the 12 3 lattice, together with the curves obtained using the parameters of fit II. The existence of a critical point separating a chirally symmetric phase from a phase where chiral symmetry is spontaneously broken can be easily spotted from the figures according to the above criterion. The corresponding lattice data can be fitted to the equation of state, yielding stable results and reasonably small χ 2 . The detailed results from the fits are reported in Tabs. 6 and 7. Fits I and II agree within errors for both values of N f on the 12 3 lattice, justifying the use of the fit with the reduced number of parameters. From a practical point of view, such a justification becomes very important since a smaller number of variables is more easily fitted from the relatively small datasets that one can produce when dealing with dynamical fermions. In the following the results quoted will always be the ones from fit II, unless specified. The critical coupling g c does not seem to be unduly affected by finite-size effects, as one can see by comparing the results in Tab. 6, 7. The values extracted from the 12 3 data for N f = 2, 4 can thus be used to check the SD prediction of equation (1.9): Comparing with our result c (N f = 4) = 0.342 ± 0.015 (3.14) we can conclude that the above prediction is more than three standard deviations away from the lattice data. This seems to favour the scenario where a critical number of flavors N f c , exists. For N f > N f c the theory would then be chirally symmetric for any value of the coupling. This conjecture is also supported by our data for N f = 6, shown in Fig. 6, where it is very difficult to identify a critical point corresponding to a symmetry breaking phase transition. Qualitatively as g increases the trajectories seem to be accumulating around a line which if continued would intercept the horizontal axis. Quantitatively, neither fit I nor fit II provides a satisfactory description of the data. It is interesting to note that, in all the cases where the fit converged, the critical exponent δ is at least two standard deviations away from the mean field value δ = 3. Finally, let us remark that finite-size effects influence the value of the critical exponent δ: the next sub-section is dedicated to the inclusion of these effects in our equation of state, so that the exponents extracted by our fitting procedure are those corresponding to the thermodynamic limit.

Finite size scaling and the equation of state
Finite size effects can be used to extract further information about the nature of the critical point. In the usual approach to finite size scaling, the RGE are assumed to be completely insensitive to finite size effects, since renormalisation deals with the short distance properties of the theory; they are therefore unchanged in a finite geometry. Nonetheless the finiteness of the volume affects the solutions of RGE, which now depend on one more dimensionful parameter. This is equivalent to considering the inverse linear size of the lattice, 1/L, as an additional relevant scaling field with eigenvalue 1 close to a fixed point at L → ∞. By repeating the arguments explained above, we obtain an equation expressing the external field as a universal function of two rescaled variables: which can again be expanded to yield a new equation of state, m = B χχ δ + A(t + CL −1/ν ) χχ δ−1/β + higher order terms, (3.16) to which the data can be fitted, in principle allowing the extraction of three critical exponents (δ, β and ν) and the critical coupling, taking into account finite size effects. The parameters obtained from this fit are the thermodynamic (infinite volume) limit ones. Hence, assuming the existence of a UV fixed point, the fit to this new equation of state provides a method to determine the thermodynamic limit of the critical exponents and coupling without further hypothesis. Once again, it is important to stress that Eq. (3.16) is based on a Taylor expansion and is expected to describe the behaviour of the data only in a restricted region around the fixed point. As in the previous section, the set of data included in the fit is selected by maximising the χ 2 /d.o.f. In order to reduce the number of parameters in the fit, δ −1/β is set to 1 again and ν is assumed to obey the hyperscaling relation (3.9), as it should be in the vicinity of a RG fixed point so that ν = (δ + 1)/3(δ −1). The number of parameters is then reduced from seven to five. Data from lattice sizes ranging from 8 to 16 have been included in the fit, whose results are reported in Tab. 8. It is clear from the values of the critical exponents that the theory is not mean-field and that the data around the critical point seem to match well the expected behaviour of a UV fixed point.
The data for the chiral condensate vs. 1/g 2 are reported in Fig. 7 and 8 for the 12 3 and 16 3 lattices respectively, together with the curves obtained from the fit. The dashed line represents the chiral condensate for m = 0 in the thermodynamic limit.
To conclude, in this section we have shown that the chiral condensate data of the N f = 2 and N f = 4 models support the hypothesis of a power-law equation of state in the critical region. The fits obtained were also consistent with the constraint δ − 1/β = 1. In passing it is worth remarking that it is not possible to obtain a good fit to our data with a mean-field equation of state with logarithmic corrections, as used for QED 4 in [35]. The critical exponents extracted for the two models are distinct, in neither case resembling those of mean field theory, suggesting that the two models have continuum limits described by distinct field theories. Of course, this in itself is not direct evidence for an interacting continuum limit at the fixed point; to establish that would require the study of lines of constant physics in the space of bare parameters g and m, and measurement of n-point functions in the critical region. No similar fit was found for N f = 6, suggesting that in this case there is no transition, and that N f c < 6 for the Thirring model.

Susceptibilities and Spectroscopy
In this section we review results of susceptibility and spectroscopy measurements, principally on the N f = 2 model studied on a 16 3 lattice with m = 0.01 (ie. unless stated otherwise), but also showing some results from a smaller lattice, and for N f = 4, for comparison. First we consider susceptibilities, generically denoted χ, which as discussed in Sec. 2 are equivalent to bound state propagators evaluated at zero three-momentum. They take large values whenever there are long-range correlations present. There are two reasons why χ is a useful quantity to study. Firstly, it can be separated into singlet χ s and non-singlet χ ns components; χ s , formed from disconnected valence fermion lines, is thus a useful measure of the importance of sea fermion loops in the vacuum. Since there are no dynamical bosons in the original model, any long-range interactions must be mediated by fermion "bubbles". The relative strength of χ s versus χ ns is the simplest probe of this effect. Secondly, if we assume that fermion -anti-fermion composites resemble canonical bosons, then the relation holds, where M is the bound state mass and Z is a renormalisation constant associated with the projection of the composite operator onto the propagating state. Hence a large susceptibility signals a small bound state mass. However, in a non-confining theory the assumptions leading to (4.1) may not be justified, eg. when there is a contribution from the fermion -anti-fermion continuum as in the scalar channel of the three dimensional GN model [3]. In this case extraction of the bound-state mass is model dependent, whereas the susceptibility is still unambiguously defined. First we consider the scalar channel, plotting in Fig. 9 χ ls , χ lns , and χ l versus 1/g 2 (see Eq. (2.25)). In the symmetric phase the signal is dominated by χ lns , which peaks at 1/g 2 ∼ 2 − 2.2, indicating a possible light scalar state, before falling away in the broken phase. The singlet contribution is almost absent in the symmetric phase, but rises across the transition to dominate the signal by 1/g 2 = 1.6. The total longitudinal susceptibility is approximately constant in the broken phase.
Next we compare longitudinal and transverse susceptibilities (2.25,2.26), showing data for N f = 2, m = 0.01 on a 16 3 lattice in Fig.10 and for N f = 4, m = 0.02 on a 12 3 lattice in Fig. 11. The main features are that χ l ≃ χ t deep in the symmetric phase, but that χ t rises steeply in the broken phase corresponding to the appearance of a massless Goldstone pion. Note that the non-singlet contribution χ tns saturates the full Ward identity result within errors; this is as expected since singlet -nonsinglet degeneracy could only be lifted by an axial anomaly, which is absent both for staggered fermions and for three dimensions.
There is an interesting relation between the ratio R of longitudinal to transverse susceptibilities and the critical scaling behaviour of the model, first noted in [22] and subsequently developed in [35]. Using the action (2.1) and the Ward identity (2.26), we have Therefore R can be calculated using the proposed equation of state (3.11): In general R depends on m, but its limiting forms in the chiral limit are easily evaluated. In the broken phase, χχ remains non-vanishing in the chiral limit, whereas in the symmetric phase χχ δ vanishes faster than m; hence The near-degeneracy of χ l and χ t in the symmetric phase from 1/g 2 = 2.4 onwards seen in Fig. 10, and from 1/g 2 = 1.2 onwards in Fig. 11, therefore provides additional support for our preferred value of δ − 1/β = 1 in the equation of state (3.11). Exactly at the critical coupling, however, R becomes independent of m, and we have In Tab. 9 and Fig. 12 we plot R for various values of m and g near the critical point, in particular for four values of m evaluated at our best estimate of the critical coupling 1/g 2 = 1.92. It is worth noting that the ratio χ ls /χ lns exceeds 40% in this region, showing the importance of including the disconnected contributions in this particular calculation. However, we also find that R decreases significantly with m, at first sight implying that the model is in the symmetric phase at this coupling; the variation of R with m is less pronounced at 1/g 2 = 1.90, which would thus appear to be closer to the critical coupling. This apparent contradiction is resolved if we take the finite-size scaling form (3.16) seriously. Using the fit parameters from Tab. 8 the coefficient multiplying the χχ δ−1/β term actually vanishes for 1/g 2 ≃ 1.88. We should therefore expect to observe a constant R at slightly stronger couplings than the critical one on a finite system. Although we have insufficient data to verify this directly, it is clear that this kind of analysis has the potential to provide an independent estimate of g 2 c and δ. An optimistic reader would consider the near-agreement found here as a signal of the fact that the RG structure underlying the equation of state approach is indeed realised in the Thirring model. In Fig. 13 we plot the vector susceptibilities χ vs and χ vns (2.27). The data is quite noisy; the trends are that χ vns appears roughly constant in the symmetric phase and fall away in the broken phase, and that the ratio χ vs /χ vns lies between 0.2 and 0.3 throughout the region.
Next we discuss spectroscopy, which involves the calculation of the timesliced correlators C(t) defined in Sec. 2. We follow the usual procedure of extracting the mass of the lowest-lying state contributing to the two-point functions from their single exponential decay for large time separations. In performing our minimum χ 2 fits, we include correlations between the data at different timeslices. We used the bootstrap technique to estimate errors on the fitted parameters. We present results from 8 2 × 16 and 16 3 lattices.
First consider the fermion, pion and scalar channels. The fermion, since it cannot mix with any excited states, shows the clearest signal and can be fitted over the entire range by the form where µ R is the physical fermion mass and the minus sign between the forward and backward terms is due to our choice of antiperiodic boundary conditions in the timelike direction. All the quoted results are obtained from fitting in the range t ∈ [2, 13] with a typical χ 2 /dof ∼ 1. Varying the fit range produces negligible variations in the fitted parameters and χ 2 /dof. Both the scalar and pion channels were fitted by the form The pion exhibits single exponential decay from early timeslices and it is possible to fit the data in the range t ∈ [3,13] in almost all cases, again with a typical χ 2 /dof ∼ 1.
Similarly to the fermion case, fitting over a smaller time interval within this range produces minimal variations. The scalar channel is more problematic. Firstly the signals exhibit greater statistical noise. Secondly, as distinct from the pion channel, we can only fit the correlators in a restricted time interval, suggesting contamination of the signal by higher-mass states. On the 8 2 × 16 lattice it is not possible to reliably extract the ground-state signal for m = 0.05, while for m = 0.02, fits are possible but quite sensitive to the fitted range. In this latter case we therefore choose a conservative fit range t ∈ [4,11]; fitting outside this range leads to a much higher χ 2 /dof. The fits on the 16 3 data are more stable for all values of the fermion bare mass. In the strong coupling regime, where the scalar mass gets larger, it proves harder to extract the ground-state mass regardless of bare mass or lattice size, resulting in larger errors on the fitted masses in this phase. In order to improve on this behaviour, it would be necessary to move to larger lattices. With regard to this last point, our present results should therefore be considered as an attempt to find qualitative evidence of the phase structure of the model, rather than a precise evaluation of its spectrum. The detailed results are presented in Tabs. 10, 11, 12, 13 and 14. Examples of the fermion, pion and scalar signals are shown in Figs. 14, 15 and 16 respectively. A comparison between Tabs. 11 and 13 enables a crude assessment of finite volume effects. In the smaller volume the fermion mass µ R is shifted upwards significantly in the symmetric phase, while remaining relatively unaffected, within errors, deep in the broken phase. The pion, in contrast, is lighter in the smaller volume, the effect again being more pronounced in the symmetric phase, but still significant in the broken phase. The scalar mass seems the least affected by the finite volume.
In Fig. 17 we plot the measured masses against 1/g 2 for the 16 3 system at m = 0.01. While the pion mass stays constant, the fermion mass is small in the symmetric phase, but rises steeply across the transition to become heavier than the pion in the broken phase. Just on the symmetric side of the transition we have M π ≃ 2µ R , indicating that the pion is a weakly bound state -unfortunately by 1/g 2 = 2.4 this picture breaks down, perhaps due to finite volume effects. The scalar mass, meanwhile, is just greater but of the same order as the pion mass in the symmetric phase, but rises across the transition to become appreciably larger. The scalar and pion masses are respectively related to the longitudinal and transverse susceptibilities: The chiral Ward identity (2.26) then yields: As long as chiral symmetry is broken, this equation states that the pion mass squared is a linear function of the bare mass and hence vanishes when the breaking is dynamical (Goldstone theorem). Fig. 18 shows the behaviour of the pion mass squared as m → 0 in the broken phase for 1/g 2 = 1.8. The line is a least squares fit; the fitted intercept lies two standard deviations from the origin, which is another signal of systematic error. It is satisfying that the spectrum reflects the phase transition so sharply, moreover in a manner entirely consistent with standard expectations -there is no evidence for the unconventional spectrum in the symmetric phase, associated with an infinite order phase transition, discussed in [15]. Of course, it would be valuable to refine this initial study both by working closer to the chiral limit, and with more lattice spacings in the timelike direction to improve the mass estimates -clearly at present we are unable to make reliable estimates of binding energies, and hence are unable to make detailed comparisons with the predictions of the 1/N f expansion in the symmetric phase.
Finally we turn to the vector channels, where we have employed both local (2.33) and one-link "conserved" (2.34) interpolating operators. In either case we found that the signal had a strong alternating component, exemplified by the plot shown in Fig. 19. Therefore fits of the form (2.35) were used. Superficially, the main difference between the two different vector channels was that the local vector gave numerically larger but noisier signals. For both operators we experimented with two parameter fits, with B set equal to −A and M d equal to M a in (2.35) for t ∈ [3,13], then three parameter fits in which A and B are allowed to differ, and finally the most general four parameter fit in which A, B, M d and M a are all independent, the last two forms being fitted for t ∈ [2,14]. The results are summarised in Tabs. 15, 16, and plotted in Fig. 20.
The main trends we found are: (i) For all fits the mass of the local operator lies significantly above that of the conserved operator. This should be no surprise, since as discussed in Sec. 2, and displayed in Tab. 5, the two operators have different quantum numbers in three dimensions.
(ii) The masses are approximately constant throughout the symmetric phase; across the transition the signal becomes smaller and thus noisier, consistent with the falloff in χ vns of Fig. 13. Our mass fits are of limited use for 1/g 2 ≤ 1.8. In the conserved channel, however, there is evidence for a sharp rise in mass at 1/g 2 = 1.9. Hence there may be a light vector particle in the symmetric phase, in accordance with the 1/N f expansion, but not in the broken phase. (iii) Although the two parameter fits gave acceptable results, there is tentative evidence in the conserved channel that the four parameter fit is significantly better, with M a < M d , particularly from the results with m = 0.02 at 1/g 2 = 2.4. An examination of Tab. 5 reveals that the direct channel has the vector spin quantum number form ψγ µ ψ while the alternating channel has the "axial vector" formψγ 5 γ µ ψ. Inspection of the lattice interaction in the form (2.12) shows that there is a contact interaction between axial vector currents of equal strength to that between vector currents; therefore it is reasonable that the one-link operator projects onto both states with comparable strengths, resulting in C(t) almost vanishing on even timeslices (Fig. 19). Thus we find both light vectors and light axial vectors in the symmetric phase. We have no explanation for the light axial vector. Following the analysis of [8], we can examine the pole condition p 2 = −M 2 A in the axial vector channel to leading order in 1/N f : One readily sees that this equation has no bound state solution of the form M A = 2m − ε. It would be interesting to see whether this persists at next-to-leading order in 1/N f . To summarise, the results in the vector channel present the only unexpected result; the light state in the axial vector channel is the first hint of departure from behaviour predicted by the 1/N f expansion in the symmetric phase. In future studies at smaller bare fermion masses, it would be interesting to monitor the ratio χ vs /χ vns , to see whether the channel becomes dominated by fermion bubble chains as in the 1/N f expansion. It would also be interesting to repeat this analysis in simulations of noncompact lattice QED 3 to see whether the light axial vector is present or not -in this case since the A µ field is now a dynamical variable, it is possible that short-wavelength fluctuations will be suppressed by a factor 1/k 2 , preventing the strong coupling to the axial vector channel occurring in (2.12).

Conclusions
In conclusion we briefly summarise our main findings, and suggest further directions to explore. For N f = 2 and N f = 4 the Thirring model exhibits a phase transition at finite critical coupling g 2 c to a phase in which chiral symmetry is spontaneously broken. This contradicts the naive 1/N f expansion, but supports the predictions of non-perturbative methods such as the Schwinger-Dyson approach. In the vicinity of the critical point the model's equation of state seems well fitted by a renormalisationgroup inspired form corresponding to power-law scaling. The picture remains con-sistent for N f = 2 even once a finite volume scaling analysis, using data from larger lattices and smaller bare masses, is employed. The critical exponents we find differ significantly between N f = 2 and N f = 4, suggesting that the two models lie in different universality classes. Although it is tempting to infer the existence of nonperturbative continuum limits at these points, it would be nice to have direct evidence by measurement of a renormalised coupling constant via a 3-or 4-point correlation function.
For N f = 6 the picture is not so clear, and we did not find a fixed-point fit for the equation of state. The Fisher plot of Fig. 6 strongly suggests that in the chiral limit the model lies in the symmetric phase for all values of the coupling, and hence that N f c < 6. However, bearing in mind the practical difficulties in identifying the strong coupling limit discussed in Sec. 2, perhaps a more conservative conclusion is that the N f = 6 model seems qualitatively very different from N f = 2, 4 [33]. It might be very difficult ever to exclude absolutely a chiral condensate exponentially suppressed in N f for large N f ; however our data disagrees with the only theoretical prediction (1.9) to put forward this scenario [11]. We intend to extend the N f = 6 studies to larger lattices and smaller m to pursue this issue further.
Although we made no attempt to test the scaling form (1.8), corresponding to an infinite order phase transition, our studies of the susceptibilities and spectrum prejudice us against this scenario: the relative ordering of the physical fermion, scalar and pseudoscalar masses varies continuously across the transition; there appear to be light bound states in the symmetric phase, comparable with the scale set by the physical fermion mass, in contradiction to the behaviour expected near a conformal transition discussed in Refs. [15,16]. It would be interesting to extend this study to N f = 4, where g c lies closer to the strong coupling limit, to see if the spectrum shows any qualitative difference. Work in this direction is in progress.
Finally, although our admittedly exploratory spectrum analysis has yielded important information, it is clear that larger lattices in the timelike direction will be needed before accurate measurements of masses and binding energies, or indeed departures from the exponential form of the decay corresponding to a spectral function more complicated than simple isolated poles, can be accomplished. Accurate mass measurements would be desirable in order to compare with predictions from the 1/N f expansion, to check whether the latter approach has any validity near the transition even in the symmetric phase. There is a first hint of its breaking down in the apparent existence of a light axial vector state here; however it is clear much more work is needed both numerically and analytically.

Acknowledgements
LDD is supported by an EU HCM Institutional Fellowship, SJH by a PPARC Advanced Fellowship. Some of the computing work was performed using resources made available under PPARC research grant GR/J67475, and those of the UKQCD collaboration under GR/K41663, GR/K455745 and GR/L29927. We have enjoyed discussing aspects of this work with Roger Horsley, Yoonbai Kim, Kei-Ichi Kondo and Nick Mavromatos.     Table 5: Spin/flavor assignments of the fermion bilinears used in spectroscopy direct alternating pion                       Figure 12: Susceptibility ratio R vs. m for 1/g 2 = 1.9 (circles), 1.92 (squares) and 2.0 (diamonds). The horizontal band marks the value of 1/δ from the fit (3.16).       Figure 19: Conserved vector correlator, plus 2 parameter fit, at 1/g 2 = 2.4.