The Chiral Critical Point in 3-Flavour QCD

We determine the second order endpoint of the line of first order phase transitions, which occur in the light quark mass regime of 3-flavour QCD at finite temperature, and analyze universal properties of this chiral critical point. A detailed analysis of Binder cumulants and the joint probability distributions of energy like and ordering-field like observables confirms that the chiral critical point belongs to the universality class of the three dimensional Ising model. From a calculation with improved gauge and staggered fermion actions we estimate that the transition is first order for pseudo-scalar meson masses less than about 200 MeV.


Introduction
In the chiral limit the order of the QCD phase transition depends on the number of quark flavours that become massless. It has been conjectured [1] and verified in numerical calculations [2,3,4] that this phase transition is first order for QCD with n f ≥ 3 massless quark flavours. The strength of this first order transition weakens in the presence of non-vanishing quark masses and eventually ends in a second order endpoint. The analysis of effective models constructed in terms of the chiral order parameter suggests that this endpoint, the chiral critical point, belongs to the universality class of the three dimensional Ising model [5]. In fact, the entire surface of second order phase transitions which in the n f -dimensional space of quark masses separates the regime of first order transitions at small quark masses from the crossover region is expected to belong to this universality class.
In the case of degenerate quark masses the lattice formulation of QCD depends on two bare couplings, the bare quark mass m and the gauge coupling β ≡ 6/g 2 . The bare quark mass acts like an external symmetry breaking field and thus leads to an explicit breaking of the SU L (n f ) × SU R (n f ) chiral symmetry. For n f ≥ 3 there is a first order phase transition in the chiral limit (m ≡ 0), which will continue to persist for small but non-zero values of m up to a critical value,m, of the quark mass. It is expected that the universal properties of this chiral critical point are controlled by a global Z(2) symmetry, which however is not an obvious global symmetry of the QCD Lagrangian. It rather is the relevant symmetry of the effective Hamiltonian which controls the critical behaviour at this point. Although the transition is second order at the chiral critical point, it is obvious that neither the chiral condensate nor the Polyakov loop will be an order parameter for the spontaneous Z(2) symmetry breaking at this critical point. In fact, we do not know a priori what are the relevant observables representing the energy-like and ordering-field like operators of the effective Hamiltonian and the corresponding scaling fields (couplings) that define the energy-like and ordering-field like directions in the coupling parameter space; it is part of the problem of analyzing the universal behaviour at this endpoint that one has to identify these operators and couplings appropriately.
The problem of determining the critical properties at the endpoint of a line of first order phase transitions is well known from studies of other statistical and field theoretic models. The construction of appropriate scaling fields has been discussed in detail in the case of the liquid-gas phase transition [6]. The concepts developed in this context have recently also been used to locate and explore the properties of the critical endpoint of the electro-weak phase transition [7,8] as well as the critical point in the ferromagnetic, three dimensional 3-state Potts model [9]. The latter problem, of course, is closely related to the analysis of the endpoint of the line of first order phase transitions that occur in the heavy quark mass limit of QCD [10].
We will focus here on a discussion of the light quark mass regime of QCD. In particular we concentrate on an analysis of universal properties of the chiral critical point in 3-flavour QCD. Although we will eventually also be interested in an accurate determination of the location of this endpoint in terms of physical mass values, e.g. quark or hadron masses, our main emphasis here will be on the analysis of universal properties. We thus will mainly use standard (unimproved) Wilson gauge and staggered fermion actions for our analysis. Lattice calculations will be performed on lattices with temporal extent N τ = 4 and various spatial lattice sizes to use well known finite size scaling methods at the critical point. A first estimate for the location of the critical endpoint on lattices with temporal extent N τ = 4 has been given in [11]. We confirm this estimate here and will, in addition, present a detailed analysis of the universal properties of the transition at the chiral critical point. Using improved gauge and staggered fermion actions we then will give a first estimate of the pseudo-scalar meson mass at the chiral critical point. This paper is organized as follows. In the next section we introduce the basic variables relevant for the description of critical behaviour in the vicinity of the chiral critical point. In section 3 we locate the chiral critical point and determine its universality class using Binder cumulants. In section 4 we construct the order parameter and energy-like operator at the chiral critical point. Section 5 is devoted to the determination of the critical pseudo-scalar meson mass at the critical point. Finally we give our conclusions in Section 6.

Scaling fields at the chiral critical point
The thermodynamics of QCD with n f quark flavours is described in terms of the partition function which is defined on a 4-dimensional lattice of size N 3 σ × N τ . Here β = 6/g 2 and m denote the gauge coupling and bare quark mass, respectively. The Euclidean action, S, is given in terms of the pure gauge action S G and the fermion matrix Q F . Throughout this paper we will use staggered fermion actions and standard Wilson or tree level improved gauge actions. The discretized QCD action then reads, The generic phase diagram of QCD in the small quark mass region is shown in Fig. 1 in the plane of bare couplings appearing in the QCD partition function. In the vicinity of the second order endpoint, (β c (m),m), the dynamics and the universal critical behaviour is controlled by an effective Hamiltonian, which can be expressed in terms of two operators E, M, i.e. the energy-like and ordering-field like operators that couple to two relevant scaling fields τ and ξ, Under renormalization group transformations the multiplicative rescaling of the couplings ξ and τ is controlled by the two relevant eigenvalues that characterize the universal critical behaviour in the vicinity of the second order critical point. The singular part of the free energy density thus scales like where the dimensionless scale factor b ≡ LT = N σ /N τ gives the spatial extent of the lattices in units of the inverse temperature. Susceptibilities constructed from E and M will show the standard finite size scaling behaviour Here V = N 3 σ denotes the 3-dimensional spatial volume and δX ≡ X − X . Note that E and M as well as S G andψψ are defined to be extensive quantities.
For finite values of the quark mass all global symmetries of the QCD Lagrangian are explicitly broken. As the symmetry of H ef f that characterizes the critical behaviour at the chiral critical point is not shared in any obvious way by the QCD Lagrangian we also may expect that in the vicinity of the chiral critical point the operators appearing in the QCD Lagrangian are mixtures of the energy-like (E) and ordering-field like (M) operators. Similarly the couplings appearing in the QCD Lagrangian are linear combinations of the scaling fields as we have indicated in Fig. 1.
In the vicinity of the critical point one may use a linear ansatz for the couplings as well as for E and M which are constructed in terms of operators appearing in the original QCD Lagrangian, Hereψψ denotes the chiral condensate evaluated on a given gauge field configuration {U}. In terms of the staggered fermion matrix, Q F (m), this is given byψψ As the operators of the QCD Lagrangian, e.g. S G andψψ or related observables like the Polyakov loop expectation value, are mixtures of E and M, the corresponding susceptibilities will all receive contributions from fluctuations of E as well as M. Asymptotically therefore all of them will show identical finite size scaling behaviour which will be dominated by the larger of the two exponents α/ν and γ/ν, respectively. For the symmetry groups of interest in the QCD context, e.g. the symmetry of three dimensional Z(2) or O(N) spin models, this will be γ/ν. A finite size scaling analysis of susceptibilities constructed from the basic operators of the QCD Lagrangian thus will give access only to the ratio γ/ν, which unfortunately is quite similar for all the above mentioned symmetry groups and thus is not a good indicator for the universality class controlling the critical behaviour in the vicinity of the chiral critical point.
The situation is different for cumulants constructed from linear combinations of ψψ and S G , From Eq. 4 it follows that for arbitrary values of x the cumulants are renormalization group invariants which in the infinite volume limit take on a universal value at the critical point (τ, ξ) ≡ (0, 0). For all values of x different from 1/r the cumulants behave asymptotically like the Binder cumulant for the order parameter; cumulants calculated on different size lattices for different quark masses will intersect at some value of the quark mass. In the infinite volume limit these intersection points will converge to a universal value which is characteristic for the universality class of the underlying effective Hamiltonian and, in fact, is quite different for the classes of three dimensional Z(2) and O(N) symmetric spin models; e.g.  (3) for O(4) [14]. The cumulants B 4 (x) thus seem to be appropriate observables to locate the chiral critical point as well as to determine its universality class without knowing in detail the correct scaling fields.

Locating the chiral critical point
The determination of the chiral critical point proceeds in two steps. First of all, we determine for fixed values of the quark mass pseudo-critical couplings, β pc (m), on finite lattices. These are defined as the position of maxima in susceptibilities of ψψ, the gauge action S G and the Polyakov loop L [15]. We then make use of the finite size scaling properties of Binder cumulants B 4 (x) evaluated at β pc (m). When analyzed as function of the bare quark masses the cumulants calculated on lattices with spatial extent L 1 and L 2 will intersect at a mass m L 1 ,L 2 . For (L 1 , L 2 ) → (∞, ∞) these intersection points will converge to the chiral critical point.
From previous studies with standard Wilson gauge and staggered fermion actions one knows that the endpoint in 3-flavour QCD is located close to m = 0.035 [11]. As the universal properties of the endpoint are not expected to be influenced by lattice cut-off effects we took advantage of this knowledge and performed our detailed scaling analysis with unimproved actions. We have performed calculations on lattices of size N 3 σ × 4, with N σ = 8, 12 and 16. We have used four values of the quark mass in the interval m ∈ [0.03, 0.04] and for each of these masses we calculated thermodynamic observables for 3 to 4 different values of the gauge coupling β. In general we collected for each pair of couplings (1 − 3) · 10 4 configurations generated with the hybrid-R algorithm a . Interpolations between results from different β-values have been performed using the Ferrenberg-Swendsen multi-histogram technique [16].  In Tab. 1 we summarize our results for the peak heights of the three different susceptibilities which we have analyzed. We generally find that the positions of these peaks coincide within statistical errors. Thus in Tab. 2 only the pseudocritical couplings extracted from the location of the peak in (δψψ) 2 are given. Also given in this table is the Binder cumulant B 4 (x) calculated at the pseudocritical couplings for two values of x. We have checked that the cumulants indeed attain their minimum at β pc (m). The case x = 0 corresponds to the cumulant of the chiral condensate alone, whereas x = 0.43 corresponds to our best estimate for the mixing parameter s, whose determination we are going to discuss in the next section.
The cumulant of the chiral condensate, B 4 (0), is shown in Fig. 2. We note that the cumulants calculated on different size lattices intersect at a quark mass close to m = 0.035. The value of B 4 (0) at the intersection point is compatible with the universal value of the Binder cumulant for the 3-dimensional Ising model. It is obvious from the data given in Tab. 2 that the situation is very similar for x = 0.43. We have fitted the cumulants on a given lattice size using a linear ansatz in the quark mass, B 4 (x) = a 0 + a 1 m. From this we find for the intersection point and the cumulant,    In addition to the Z(2) value for the Binder cumulant we show in Fig. 2 also the result for 3-dimensional O(2) symmetric spin models. As the Binder cumulant depends quite sensitively on the underlying symmetry at the critical point the result found by us for the crossing point strongly suggests that the chiral critical point indeed belongs to the universality class of the 3-dimensional Ising model. The ratio of critical exponents γ/ν on the other hand is not sensitive to the universality class. As expected one finds, however, that all three susceptibilities show identical finite size scaling behaviour; for the quark mass closest to the critical point, m = 0.035, the ratio of susceptibilities calculated on lattices of size N σ = 12 and 16 takes on the value 1.69(23), 1.66(22) and 1.66(24) for chiral, Polyakov-loop and action susceptibilities, respectively. This corresponds to a ratio of critical exponents γ/ν = 1.8(5) which is consistent with the 3-dimensional Ising as well as O(2) and O(4) values (γ/ν ≃ 1.96).
As discussed above the intersection point of Binder cumulants constructed from M(x) according to Eq. 8 will be independent of the choice of x in the infinite volume limit. In fact, from Eq. 4 it follows that for x = s the singular part of the free energy will lead to a unique crossing point, independent of the lattice size, i.e. finite volume effects are minimized if we manage to select the correct ordering-field like operator for constructing the Binder cumulant. For all other choices a volume dependence of the intersection points calculated from Binder cumulants on different size lattices will result from the scaling behaviour of f s (τ, ξ) and they will converge to the universal value only in the infinite volume limit. It follows from Eq. 4 that the  volume dependence is quadratic in ∆ = x − s, i.e. in deviations from the optimal choice of x. This behaviour is confirmed through the analysis of the x-dependence of the intersection point of the Binder cumulant calculated on lattices of two different sizes. In Fig. 3 we show the intersection points of Binder cumulants calculated on lattices of size 8 3 × 4 and 16 3 × 4. As can be seen the volume dependence is small for a wide range of x values and the extremum, which corresponds to the optimal choice x = s, is closest to the 3-d Ising value. In fact, determining the x-value where the extremum is reached provides a finite volume estimate for the mixing parameter s. From a jackknife analysis one finds for the location of the minimum s min = 0.430(23).

The mixing parameters r and s
Having located the chiral critical point we now can determine the mixing parameters r and s. This will allow us to construct the energy-and ordering-field like operators and obtain further evidence that the chiral critical point belongs to the universality class of the 3-dimensional Ising model. vicinity of the critical point. We thus demand or equivalently Using these relations together with Eqs. 6 and 7 we obtain two conditions for the mixing parameters r and s in terms of the parameter B and expectation values of S G andψψ, The parameter B controls the mixing of the gauge coupling and bare quark masses needed to define lines of constant ξ. In the symmetry broken phase one line of constant ξ is known to us; the line of first order phase transitions defines the zero external field line (ξ = 0) of the effective Hamiltonian. We thus can extract B from the quark mass dependence of the pseudo-critical couplings. Using Eq. 6 one obtains Knowing B we also know r and can construct s using Eq. 12.
In order to determine the mixing parameter r from Eq. 13 we have to approximate the derivative by finite differences as we can perform calculations of β c (m) only at a discrete set of quark mass values. A first estimate may be given using our data on the largest lattice (16 3 × 4). The slope of β pc (m) defines r −1 . We estimate this from a straight line fit to the values given in Tab. 2 which gives, r = 0.51 (2). The mixing parameter r is large and definitely non-zero; as expected, a mixture of two operators is needed to construct an energy like observable in which the otherwise dominant ordering-field like contributions to S G andψψ cancel. The importance of choosing the correct mixing parameter r becomes apparent from an analysis of joint probability distributions for δE and δM. These are little affected by changes of s but they strongly depend on r. In Fig. 4 we show contour plots for the joint probability distributions and various values of r and s. For r = s = 0 these are just the fluctuations in δS G and δψψ which are strongly correlated. Only for larger values of r fluctuations in the order parameter are independent of changes in the energylike observable. In fact, we have used this criterion to improve our determination of r over our previous estimate based on the slope of β c (m) atm. We demand that fluctuations in δM vanish for any fixed value of δE. This generalizes Eq. 11 and maximizes the Z(2) symmetry of the contour plots shown in Fig.4 around the δM = 0 axis. As an estimate for the mixing parameter r one obtains in this way r = 0.550 (7). A determination of the parameter s from Eq. 12 then yields s = 0.41(51). Although this value has large errors it is consistent with the result found from the extremum of intersection point of the Binder cumulant shown in Fig. 3. As a best estimate of the mixing parameters we therefore obtain r = 0.550 ± 0.007 , s = 0.430 ± 0.023 . (14) 5 Physical scale at the chiral critical point   (1), and using estimates for the critical temperature in 2 and 3-flavour QCD [17] we thus estimate for the pseudo-scalar meson mass at the chiral critical point m ps ≃ 290 MeV.
The entire analysis of the chiral critical point discussed so far has been performed with unimproved gauge and fermion actions on rather coarse lattices. Improved actions are not expected to modify the results on the universal properties of the chiral critical point, which have been presented above. They may, however, well influence the quantitative determination of the chiral critical point. It has been found in studies of the first order deconfinement transition occurring in the pure gauge sector that the gap in physical observables like the latent heat or surface tension is cut-off dependent. Improved actions generally lead to smaller gaps and a reduced cut-off dependence of these observables [18]. One thus may expect that also in the region of first order chiral transitions the gap in the chiral condensate gets reduced when calculated with improved actions. This will shift the critical point to smaller values of the pseudo-scalar meson mass. We therefore have investigated the chiral critical point also in calculations with improved gauge and staggered fermion actions (p4-action), which we have used previously for studies of the thermodynamics of two and three flavour QCD [17,19]. In these calculations no evidence for a first order transition has been found down to bare quark masses m = 0.01 [19]. We now have extended these calculations to smaller quark masses and made use of the universal properties of the Binder cumulants discussed above. One can locate the chiral critical point quite accurately through a calculation of Binder cumulants on finite lattices and there is no need to accurately determine the correct order parameter for such an analysis; i.e. the analysis can be performed with Binder cumulants constructed from the chiral condensate. We have performed calculations on a lattice of size 12 3 × 4 with a bare quark mass m = 0.005 and on a lattice of size 16 3 × 4 with m = 0.01. Again a Ferrenberg-Swendsen reweighting is used to determine the pseudo-critical couplings and Binder cumulants at these couplings. The results are summarized in Tab. 3.
From the results given in Tab. 3 it is obvious that the smaller quark mass leads to a Binder cumulant below the value for the 3-dimensional Ising model and thus is in the first order region of the 3-flavour phase diagram. This also is confirmed by the large value of the chiral susceptibility. For the larger quark mass, on the other hand, the susceptibility stays small and the Binder cumulant is significantly above the Ising value. This quark mass still lies in the crossover region. We thus have obtained an upper and lower limit for the chiral critical point, suggesting a critical bare quark mass ofm = 0.0075(25). Extrapolating the meson masses calculated in [17] to this value of the bare quark mass, we estimate for pseudo-scalar meson mass at the chiral critical point m ps ≃ 192(25) MeV. The critical mass thus is considerably smaller than estimated from calculations with unimproved gauge and fermion actions.

Conclusions
Through an analysis of Binder cumulants we have verified that the chiral critical point in three flavour QCD belongs to the universality class of the three dimensional Ising model. The analysis of joint probability distributions provides a powerful tool to construct the order parameter at the chiral critical point as well as the energy like scaling field. Although the chiral condensate itself is not the order parameter at this critical point, we explicitly have verified that Binder cumulants constructed from it are little influenced by finite volume effects and are good observables to locate the critical point as well as the universality class.
Having determined the universality class of the chiral critical point for three degenerate quark masses one could use this information to determine the critical parameters also for non-degenerate quarks from calculations of Binder cumulants on finite lattices. The quality of the straight line fits shown in Fig. 2 suggests that a first order Taylor expansion of B 4 (0) in terms of degenerate up/down quark masses m u,d and a strange quark mass m s around the three flavour critical pointm = 0.033 might be possible. As the chiral critical line corresponds to those sets of quark masses where the Binder cumulant attains the 3-d Ising value one has to determine the line on which B 4 (0) stays constant in the (m u,d , m s )-plane. To leading order this line is given by For instance, for m u,d = 0.025 this estimate suggests that the region of first order chiral transitions ends at a critical value of the strange quark mass m s ≃ 0.049. This is consistent with the result of [4] where no sign for a first order transition has been found for (m u,d , m s ) = (0.025, 0.1) and not in contradiction with [11] who reported on a "weak first order like behaviour" at (m u,d , m s ) = (0.025, 0.05) as opposed to two state signals at lower m s values.
The chiral critical point has been determined in calculations with unimproved as well as improved gauge and staggered fermion actions on lattices with temporal extent N τ = 4. The physical scale extracted from calculations of the pseudo-scalar meson mass at these endpoints is quite different in both cases. This indicates that cut-off effects are still significant and calculations closer to the continuum limit are definitely needed to fix a physical scale for the location of the chiral critical point. The present indication is, however, that improved actions lead to smaller values for the critical pseudo-scalar meson mass. This makes it increasingly unlikely that the transition in the physically realized case of two light and a heavier strange quark lies in the first order region of the QCD phase diagram.