CERN-PH-TH/2006-081

MS-TP-06-03

{centering}

The chiral critical line of QCD

at
zero and non-zero baryon density

Philippe de Forcrand and Owe Philipsen

Institut für Theoretische Physik,
ETH Zürich,
CH-8093 Zürich,
Switzerland

CERN, Physics Department, TH Unit, CH-1211 Geneva 23,
Switzerland

Institut für Theoretische Physik, Westfälische Wilhelms-Universität Münster, Germany

We present numerical results for the location of the chiral critical line at finite temperature and zero and non-zero baryon density for QCD with flavours of staggered fermions on lattices with temporal extent . For degenerate quark masses, we compare our results obtained with the exact RHMC algorithm with earlier, inexact R-algorithm results and find a reduction of 25% in the critical quark mass, for which the first order phase transition changes to a smooth crossover. Extending our analysis to non-degenerate quark masses, we map out the chiral critical line up to the neighbourhood of the physical point, which we confirm to be in the crossover region. Our data are consistent with a tricritical point at 500) MeV.

We also investigate the shift of the critical line with finite baryon density, by simulating with an imaginary chemical potential for which there is no sign problem. We observe this shift to be very small or, conversely, the critical endpoint to be extremely quark mass sensitive. Moreover, the sign of this shift is opposite to standard expectations. If confirmed on a finer lattice, it implies the absence of a critical endpoint for physical QCD at small chemical potential, or another revision of the QCD phase diagram. We critically examine earlier lattice determinations of the QCD critical point, and find them to be in no contradiction with our conclusion. Hence we argue that finer lattices are required to settle even the qualitative features of the QCD phase diagram.

## 1 Introduction

Based on the property of asymptotic freedom, a fundamental prediction of QCD with three flavours of quarks is the transition from the familiar hadronic physics at low temperatures to a regime of “deconfined” quark gluon plasma at high temperatures. Whether this transition is characterised by singular behaviour of the partition function corresponding to a first or second order phase transition, or merely represents a smooth and analytic crossover between different dynamical regimes, depends crucially on the choice of the quark masses and the net baryon density specified by its chemical potential, . In the following we shall assume the light quarks to be degenerate, , and vary independently. We couple the quark chemical potential to the light quarks only, except for the degenerate case , where all quarks are coupled. The parameter space of the theory considered here is thus four-dimensional, .

The first task in determining the phase diagram in this parameter space consists of finding the (pseudo-)critical temperature , defined e.g. by the peak of some susceptibility, which represents the boundary between the hadronic and plasma regimes. The independent variables then span a 3d parameter space with regions of first order phase transitions and analytic crossover separated by surfaces of second order phase transitions. In order to identify the order of phase transitions, and the location of the critical surfaces in particular, finite size scaling analyses are necessary.

Let us first discuss the situation for , shown schematically in Fig. 1 (for early references, see, e.g. [1]). Gauge invariant, local order parameters characterising the transition only exist in the extreme cases of zero or infinite quark masses, namely the chiral condensate and the Polyakov loop, respectively. These limiting theories thus must feature singular phase transitions, and one may write down effective theories of the Ginzburg-Landau type for the order parameters [2]. It is numerically well-established that the phase transition is first order in the quenched [4] limit, and there is strong numerical evidence for first order in the chiral [3] limit. Since first order phase transitions are robust against small variations of the parameters of the theory, the first order regions must extend by a finite amount into the quark mass plane. On the other hand, simulations have revealed smooth crossover behaviour for intermediate quark masses, which implies second order boundary lines between the first order and crossover regions.

In the case of heavy dynamical quarks, the relevant symmetry is the center symmetry, and the weakening,
by the dynamical quarks, of the first-order transition which occurs in the Yang-Mills theory, is understood
qualitatively [5] and to a large extent quantitatively. Simulations have determined the second-order line to
correspond to a meson mass of about 2 GeV [6], and have confirmed the expectation that the universality class is that of the 3d Ising model.
In the case of light quarks, numerical simulations are more difficult, and
very little is known quantitatively about the location of the second order boundary line.
The only point computed to some accuracy with standard staggered fermions is the chiral critical point
^{1}^{1}1The superscript ”” here and in the following refers to ”critical”, not to charm.
on the diagonal [7, 8, 9], which also was determined to belong to the 3d Ising universality class [7].

While the statement about the universality class concerns infrared physics and thus is stable against cut-off effects, the location of the critical point in the physical mass plane turns out to be very strongly affected. To date calculations have been performed on lattices with 4 time-slices only (, implying a lattice spacing fm), but simulations with improved actions indicate values for , and the associated pion mass, which are considerably smaller than the standard action result [7]. Moreover, all these simulations used the so-called R-algorithm [10], which has stepsize errors and therefore gives only approximate results in the absence of a careful extrapolation to zero stepsize. In any case, all current results are consistent with the physical point being in the crossover regime.

In the presence of a chemical potential the second order boundary lines turn into surfaces, as indicated in Fig. 2. The qualitative features of the phase diagram now depend crucially on the curvature at , . The common expectation is that this curvature is positive. Hence the physical point, once the chemical potential is increased, will be closer to the critical line, and intersect it for a critical chemical potential . For values larger than a first order phase transition is expected. Clearly, this is not the case for negative curvature of the critical surface.

In this work we present a comprehensive numerical study mapping out the chiral critical line in simulations of the standard staggered action on several lattices with . Upon repeating the computation for the chiral critical point with the rational hybrid Monte Carlo (RHMC) algorithm [11], which is free of finite step size errors, we find that the bare quark mass is reduced by 25%, and the physical pion mass by 10%, compared to the accepted values determined previously using the R-algorithm. We then extend our simulations to cover a wide range of quark masses, mapping out the critical line up to the neighbourhood of the physical point. In agreement with expectations, the physical point is found to be on the crossover side of the boundary. Assuming O(4) behaviour for the chiral limit, the fit to our critical line can be extrapolated to the axis consistently with the required scaling behaviour, putting the tri-critical point in that scenario (see Fig. 1) around . However, non-O(4) behaviour is not excluded by our data. Our results should also provide a testing ground and input for analytic attempts to determine the critical line from effective theories based on universality arguments [12] (for a review, see [13]).

In a second set of simulations, we repeat the analysis for an imaginary baryon chemical potential and determine the corresponding shift of the critical line, following the strategy already used in [9]. Together with additional imaginary simulations for the case, this allows for a determination of the curvature of the critical surface at , which can be readily continued to real values of . We find this curvature to be negative, as illustrated in Fig. 2 (right). In the phase diagram this implies that the critical endpoint moves to smaller with growing quark mass, until it disappears entirely for physical quark masses. This is contrary to customary expectations, and in contradiction with the results of [14], obtained at the same lattice spacing and with the same action, but using the R-algorithm and a different numerical approach. Clearly, a careful study of systematic errors, due in particular to the very coarse lattice spacing, is needed. Still, if the physical point of QCD is indeed in the crossover region at , our finding would imply that the transition will remain an analytic crossover also for any finite MeV, placing a possible QCD critical point at much larger values of . Preliminary results to this extent have already been given in [15].

After summarising the properties of QCD at imaginary and introducing the Binder cumulant as our observable for the order of the phase transition in Secs. 2, 3, respectively, we begin our analysis in Sec. 4 with a thorough discussion of step size effects for and a comparison of results from the R- and RHMC-algorithms. The computation of the chiral critical line for and is presented in Sec. 5, which also discusses the resulting new scenario for the () phase diagram of physical QCD. An assessment of systematic uncertainties is contained in Sec. 6, along with our conclusions.

## 2 QCD at imaginary

In order to study the phase diagram Fig. 2 at finite baryon density, we employ simulations at imaginary chemical potential , where the fermion determinant is positive, followed by analytic continuation, as discussed in detail in previous work [16, 9]. To render the paper self-contained, we briefly recall some points needed in the sequel. The QCD partition function at finite baryon chemical potential is even under reflection . Moreover, it is periodic in the imaginary direction, with period for colours [17], i.e. . Because of the fermionic boundary conditions, this symmetry implies that a shift in by is exactly compensated by a -transformation, so that transitions take place between neighbouring centre sectors for all 16, 18], as in Fig. 3. Hence, the first of these transitions limits the radius of convergence for analytic continuation to the first sector for most observables. With a pseudo-critical temperature of MeV, our accessible physics range is thus MeV. . It has been numerically verified that these transitions are first order for high temperatures and a smooth crossover for low temperatures [

Within this first sector, observables can be simulated at imaginary . The results may be fitted by truncated Taylor series , whose convergence can be tested by inspection. Analytic continuation of successful fits is then trivial.

A remarkable finding of previous work inspecting Taylor series is that, within MeV, convergence is rapid and screening masses [19, 20], the pressure [21] as well as the pseudo-critical temperature [16, 9] are all well described by the leading term . This becomes plausible by noting that at finite the natural expansion parameter is rather than [19, 20]. Hence we write the pseudo-critical temperature, separating the hadronic from the plasma phase, and the critical quark mass marking the boundary between first order and crossover behaviour, as

(1) | |||||

(2) |

The choice to treat as the independent variable parametrizing the critical line in Eq. (2) reflects our practical procedure, namely to fix and then scan in , because the critical line is a steeper function of the latter. We shall determine the coefficients quantitatively and free of step size errors for the theory, and provide the sign of along the whole chiral critical line for .

## 3 Universality and the Binder cumulant

There are different ways to investigate and exhibit critical behaviour of the theory along the critical line, such as finite size scaling (FSS) of susceptibilities, of Lee-Yang zeroes or of the Binder cumulant. The Binder cumulant [22] offers various advantages for our particular study. It is defined as

(3) |

with the fluctuation of some observable around its mean value, evaluated at the pseudocritical coupling (determined by a peak in or a zero in ). In practice we study for the two quark masses . In the infinite volume limit the Binder cumulant behaves discontinuously, assuming the values 1 in a first order regime, 3 in a crossover regime and some critical value reflecting the universality class at a second order critical point. On a finite volume the discontinuities are smeared out and flattened, so that passes continuously through the critical value. The location of the critical point in parameter space will then be displaced by some finite volume correction.

In [7] the Binder cumulant was chosen as observable because its critical value for the expected Z(2) universality class is distinct from those corresponding to other symmetries like O(2), O(4) etc. In this way Z(2) scaling for the chiral critical point was clearly established in [7]. Once the universality class is ascertained, the Binder cumulant allows to approximately map out the critical line on a fixed lattice size, by scanning the parameter space for the line on which is held constant at its critical value. In practice this is best achieved by holding one quark mass fixed, and scanning in the other. In the small regime the critical line turns out to be a steep function , and thus we choose to scan in the light quark mass while keeping fixed. In the neighbourhood of a critical point can then be expanded in a Taylor series,

(4) |

with for . For there is only one mass variable in the above expression, and the mass dependence of the coefficients disappears.

For large volumes the approach to the thermodynamic limit is governed by universality. Near a critical point the correlation length diverges as , where is the distance to the critical point in the plane of temperature and magnetic field-like variables, and for the 3d Ising universality class. In practice, we first find the pseudo-critical gauge coupling for a given pair , and then compute for those parameter values. Since is tuned to always, we have . is a function of the dimensionless ratio , or equivalently . Hence one expects the universal scaling behavior

(5) |

## 4 without step size errors

The algorithm most widely used in simulations of finite temperature QCD in the staggered formulation, both standard and improved, is the R-algorithm [10], which was also employed in previous studies of the chiral critical point for [7, 8, 9]. As pointed out in [10], the “correct” usage of the R-algorithm consists of performing simulations for various choices of decreasing stepsizes, followed by an extrapolation to zero stepsize. However, in practice usually a shortcut avoiding the extrapolation is applied: for some reference value of the quark masses, choose a step size for which the step size error is smaller than the typical statistical error of the simulation. The molecular dynamics of the R-algorithm then suggests to keep the ratio of quark mass and step size constant, i.e. adjust the stepsize accordingly when the quark mass is reduced. A typical choice is , although in many cases or even have been adopted.

While this procedure has been followed successfully in the intermediate quark mass regime, it breaks down for small quark masses, where the linear relation no longer appears to hold and the step size needs to be decreased faster than proportionally. Furthermore, in a study of the QCD phase transition at finite isospin chemical potential it has recently been demonstrated that a finite step size leads to a systematic underestimate of [23]. Hence too coarse a stepsize can fake a first order transition, when the zero step size result really represents a crossover behaviour.

### 4.1 The order of the transition: R- vs. RHMC algorithm

In order to control this important source of systematic error, we have returned to our investigation of the critical point at [9], this time with the RHMC-algorithm. This algorithm has no stepsize errors and is exact. For a discussion of the algorithm and numerical test results, see [24]. Our numerical procedure to compute the Binder cumulant is as follows. For each set of fixed quark mass and chemical potential, we determine by interpolating from a range of typically 3-4 simulated -values by Ferrenberg-Swendsen reweighting [25]. For each simulation point 50k-200k RHMC trajectories have been accumulated, measuring the gauge action, the Polyakov loop and up to four powers of the chiral condensate after each trajectory. Thus, the estimate of for one mass value consists of at least 200k, and the estimate of the critical mass of at least 800k trajectories.

Fig. 4 (left) shows results of this first study, comparing measurements of from the RHMC algorithm with those obtained from the R-algorithm at various step sizes. The figure confirms the finding from [23] that decreasing the step size leads to an increase in the values of the Binder cumulant. It also constitutes a useful test of the RHMC algorithm, whose results indeed correspond to the zero step size limit of the R-algorithm. We note that for our smallest quark masses studied, , the RHMC algorithm runs over 20 times faster than the R-algorithm at the commonly applied step size. Since the latter also requires runs at several step sizes for the extrapolation, the RHMC is thus considerably more economical in producing results free of step size errors.

### 4.2 The critical quark mass at

In order to eliminate step size errors, we now proceed to repeat the calculation of the critical quark mass in the three-flavour theory by means of the Binder cumulant, this time with the RHMC algorithm. The result obtained on an lattice is shown in Fig. 4 (right). Qualitatively the behaviour is the same as previously, with growing from first order behaviour through its critical Ising value to crossover with increasing quark mass, which can be fitted to leading order in the quark mass. However, the critical Ising value is now obtained at a bare mass of , which is about 25% smaller than the value quoted by all previous work using the R-algorithm [7, 8, 9].

One may ask whether this change affects bare quantities only, while the R and RHMC algorithms probe the same physics. To study this issue, we measured the zero-temperature hadron spectrum at the parameters using RHMC, and compared with the same exercise performed in [7] at the parameters using the R-algorithm. For the pion, which is the most accurately determined, the ratio changes from 1.853(1) (R [7]) to 1.680(3) (RHMC). This reduction of 10% in the pion mass corresponds to a change of 20% in the renormalized quark mass, very near the observed 25% change in the bare quark mass. Therefore, replacing the R by the RHMC algorithm corrects a large error in the physical values of the critical parameters. The correction should be even larger for smaller masses. We conclude that for the study of the QCD phase transition in the region of physical quark masses, step size errors in the Monte Carlo algorithm can lead to a qualitatively different picture at fixed parameter values, and the use of an exact simulation algorithm is mandatory.

Our results so far have been obtained on a single spatial volume . The next task is to study the FSS behaviour and uncover possible finite volume effects. This is particularly important since large finite size corrections were reported in a recent investigation of the chiral critical point at finite density in the Taylor expansion of susceptibilities [26]. Fig. 5 shows data obtained for three lattice sizes with . According to Eq. (5) and the corresponding discussion, in the scaling region near a critical point should be described by

(6) |

with and for 3d Ising universality. We have checked for finite volume effects by fixing to the Ising value and fitting for and . With a of 0.74 per d.o.f., we obtain consistent with our result from only, and , which is consistent with the Ising exponent. We conclude that for the Binder cumulant is close to thermodynamic scaling for lattice sizes , and hence finite volume effects are under control in this calculation.

### 4.3 The pseudo-critical temperature for at finite

On the lattice, is determined from the pseudo-critical gauge coupling, which we define as the location of the peak of the plaquette susceptibility. On any finite volume it can be expanded as a double series in mass and chemical potential around the three flavour critical point ,

(7) |

Fig. 6 shows the measured values of for four different imaginary chemical potentials spanning the whole range up to the first -transition. For each value of , four quark masses in the range have been simulated.

5. | 1369(3) | 0. | 781(7) | — | 1. | 94(3) | 1. | 28 | |
---|---|---|---|---|---|---|---|---|---|

5. | 1371(3) | 0. | 759(22) | 0. | 33(32) | 1. | 95(4) | 1. | 27 |

As in our previous study, we obtain good fits retaining the leading terms only, as shown in Table 1. In particular, the term quadratic in the chemical potential is now sufficient to describe the data all the way out to , the quartic coefficient being consistent with zero. Using the two-loop beta function, this translates to a pseudo-critical temperature of

(8) |

Again, we note a shift of up to 10% in the coefficients compared to the R algorithm results [9]. Eq. (8) has to be considered with some caution, since it is well known that the two-loop beta function is rather inaccurate at our coarse lattice spacing. The effect of the non-perturbative beta function is to increase the absolute values of the coefficients , perhaps by up to a factor of 2 [27].

### 4.4 Quark mass dependence of the critical point for

In order to detemine the order of the transition, we now repeat the previous procedure to find how the critical bare quark mass changes with imaginary chemical potential. As in the case of the pseudo-critical temperature, we express this by a Taylor series Eq. (2):

(9) |

in order to be able to continue to real . Inversion of this function will then give the location of the critical point as function of the quark mass, . In practice, at a finite lattice spacing we are dealing with the expansion in lattice units,

(10) |

A crucial point is that, for fixed temporal lattice extent , the lattice spacing entering the dimensionless and is different, since depends on . The relation of the leading coefficient to its continuum counterpart is thus given by

(11) |

or in terms of from Eq. (8) and from Eq. (10),

(12) |

The coefficients are extracted from our data for obtained at imaginary , by fitting to a double expansion about the known critical point at ,

(13) |

The leading coefficients are then obtained as

(14) | |||||

(15) | |||||

For the actual analysis it is thus convenient to reparametrise the second order expansion of as

(16) | |||||

with , and fit the data via the parameters .

Our data for five different values of imaginary chemical potential are shown in Fig. 7 (left). Remarkably, there seems to be negligible influence of the chemical potential. The results of various simultaneous fits of all four curves are displayed in Table 2. All fits are good, and none of the next-to-leading terms is significantly constrained. This is corroborated by discarding all next-to-leading terms, which leads to a perfectly acceptable fit with consistent with zero, as in the last line of Table 2. Fig. 7 (right) displays the error band coming from a linear fit (Table 2, line 3). Clearly, the slope is very nearly zero.

0. | 0262(7) | 13. | 3(1.4) | -91. | 6(143.5) | -0. | 079(47) | -1. | 6(1.0) | -2. | 1(3.5) | 0. | 90 |

0. | 0263(6) | 13. | 9(0.6) | — | -0. | 075(42) | -1. | 35(0.73) | — | 0. | 82 | ||

0. | 0270(5) | 13. | 6(0.6) | — | -0. | 0024(160) | — | — | 0. | 93 | |||

0. | 0271(3) | 13. | 6(0.6) | — | — | — | — | 0. | 88 |

The final result is then obtained by employing Eq. (12) to convert to continuum units. The second factor in Eq. (12) is 1.077(2), close to 1. In the first factor, the term , which describes the variation of with real and is thus negative, reinforces the negative trend of , to yield

(17) |

Hence, we arrive at the surprising result that the first order region in the phase diagram Fig. 2 shrinks when a real chemical potential is switched on. This is contrary to the expected qualitative behaviour.

The reader will notice the large error on the coefficient in Eq. (17). It is a conservative estimate and stems entirely from the larger error on . If one were to include a -term, the previous conclusion would only be strengthened: the leading term gets more negative and a negative quartic term comes on top of it, cf. Table 2.

However surprising, our findings agree with preliminary results for the same lattice theory at finite isospin chemical potential, which indicate that there too, the transition becomes weaker as the chemical potential is turned on [28]. Finally, let us note that the same qualitative behaviour applies to the first order region in the heavy quark limit, which has recently been shown to also shrink with real chemical potential [29].

## 5 The chiral critical surface for

Having removed finite step size errors from the calculations, we proceed to map out the chiral critical line for non-degenerate quark masses. All our simulations have been performed with the RHMC algorithm. Since lattices proved to be large enough for our observable in the case of , we use that lattice size to trace out the critical line, performing another check of finite volume effects at on a lattice.

With two different quark masses in the theory, a technical question concerning the Binder cumulant arises. Obviously, can be constructed from the chiral condensate of either mass flavour. Universality guarantees that, in the infinite volume limit, either choice tends to the same universal value. However, in a finite volume there are corrections, and they are different for different operators. The corrections are minimised for that superposition of operators, which corresponds most closely to the mapping of the QCD parameters onto the scaling fields of the effective 3d Ising model. It is well known that, even for the case of three degenerate flavours, this is a superposition of the chiral condensate and the plaquette [7], as well as higher dimension fermionic and gauge condensates.

Here, we do not attempt to construct an optimised observable by mixing in gauge condensates, but simply compare the behaviour of constructed from condensates of different mass flavours, , as shown in Fig. 8. The observables constructed from the condensates of the heavy and light flavours, respectively, are observed to intersect the critical value at significantly different values of . Nevertheless, comparison of results obtained at shows that this difference is rapidly disappearing on larger volumes. Moreover, the common intersection point to which the results converge appears to be close to that obtained from on , indicating that the latter has far smaller finite volume corrections. This is not too surprising, as one would expect the scaling field corresponding to chiral symmetry to be dominated by the lightest quark flavour. Hence, in the following we will always work with constructed from light quark condensates.

### 5.1 The critical line for

By fixing and scanning in (at least 4 values), the critical light quark mass for that choice of is determined by interpolation, analogously to the three-flavour case. This is repeated for other values of , resulting in the sequence of critical points displayed in Fig. 9, left. As in the three-flavour case, every critical point appearing in this figure consists of at least 800k RHMC trajectories.

There are several features of Fig. 9 worth discussing. An interesting observation concerns the behaviour of the function in the neighbourhood of the three flavour critical point. If is constructed from gauge condensates and neglecting the change in the pseudo-critical coupling with the quark masses, a Taylor expansion around the symmetric critical mass yields for the line of constant (critical) the leading order result [7]

(18) |

i.e. the critical line should pass through the symmetrical point with slope -2. In contrast, our data extracted from exhibit a different slope, see Fig. 9 (left). This underlines again the importance of choosing an appropriate observable for finite volume computations. A Taylor expansion of as in Eq. (18) is only possible on finite volume, but expanding would yield additional non-perturbative contributions to Eq. (18). We thus conclude that Eq. (18) does not describe the critical line, not even in the immediate neighbourhood of .

Another interesting question is how the critical line continues to even smaller light quark masses. If the chiral limit of the theory exhibits O(4) universality, then the critical line hits the axis in a tricritical point at some finite strange quark mass value [2]. Whether this scenario is realized or not is an issue not yet settled (cf. the discussion and references in [15]). Among the most recent publications using staggered fermions, one favors a first order scenario for the chiral limit [30] while the other supports the O(4) scenario [31]. With our current data, we are unable to decide this question, but we can check for consistency with the O(4) scenario, which implies mean-field exponents near the tricritical point . Indeed our data support a fit to the approach to the chiral limit, as shown in Fig. 9, left, predicting the tricritical point to be at or . Note however that our lattice is very coarse ( fm), and our spatial volume becomes rather small as is reduced: for the uppermost point in Fig. 9, left, corresponding to the physical strange quark mass, only. Thus, our systematic error might be rather large. Nevertheless, we have strong indications that is significantly larger than the physical strange quark mass.

### 5.2 The chiral critical line and the physical point

The most important question regarding the critical line is, of course, its location relative to the physical point of QCD. So far, all known lattice data are “consistent” with the physical point being in the crossover region. This is also the result found by a simulation at physical quark masses in [14]. However, these results were obtained by the R-algorithm, and we have seen in the three flavour case that significant shifts in the critical quark masses can arise due to step size errors.

In order to estimate the location of our critical line in physical units, we have therefore performed zero temperature simulations with bare quark masses corresponding to two points on the critical line. One corresponds to the three flavour theory and the other to the point with roughly physical strange quark mass, as indicated by the arrows in Fig. 9 (left). The parameters of the simulations are given in Table 3, together with the measured meson masses. In both cases, the lattice size was , and about 400 configurations were analyzed.

(0.0265,0.0265) | 5.1374 | 0.420(1) | 0.420(1) | 1.383(7) | 0.304(2) | 0.304(2) |

(0.005,0.25) | 5.1857 | 0.2109(1) | 0.8915(1) | 1.398(16) | 0.151(2) | 0.638(8) |

Setting a physical scale along the critical line is a tricky problem. Neither of our simulation points matches the physical point, so that strictly speaking one cannot match to any real world observable. Doing so anyway, different observables inevitably give different values for the lattice spacing. A measurement of the force via elongated Wilson loops gives and 1.87(2) respectively, where fm is the Sommer scale. This amounts to fm in both cases. On the other hand, matching the -mass to its physical value gives a lattice spacing which is by 20% larger. Note, however, that a difference of similar magnitude has been observed in [14] on the physical point. This suggests that the greater part of this difference is due to cut-off effects rather than to the deviation from physical parameters.

It thus appears safer to avoid setting an absolute scale altogether, and instead compare the meson mass ratios from Table 3 with their values at the physical point, and . We thus conclude that our point on the critical line indeed corresponds to the physical Kaon mass and to pions lighter than physical. In other words, the physical point is on the crossover side of the critical line.

It is interesting to compare our meson masses, in
lattice units, with those of [14], which used the same strange quark
mass, and the R-algorithm at almost the same inverse coupling .
The comparison is shown in Fig. 10, where the
straight lines are fits to the data of [14] only. Good consistency
is apparent, showing that the R-algorithm step size error for the meson
masses is small, for the parameters considered in [14]
^{2}^{2}2
Note also, that there are preliminary results by Z. Fodor and S. Katz
(https://www.bnl.gov/sewm/) using finer lattices and the exact RHMC algorithm,
which confirm that the physical quark masses give a finite-temperature
crossover.. The figure also shows to be practically independent of
, thus affirming our conclusion above.

Finally, the fact that the lattice spacing varies little between our two simulation points, implies that itself does not change much as one moves along the critical line. This is in agreement with model calculations [13].

### 5.3 The critical line for and the critical surface

We have also run a second set of simulations with an imaginary chemical potential , in order to determine how the critical line shifts with baryon density. These data are shown in Fig. 9 (right), in lattice units. We observe that the shift in the critical line is very small, despite the sizeable value of the chemical potential. Within two sigma the line is consistent with its counterpart. Moreover, to the extent that there is a displacement, it shows a trend to lie to the right of the zero density line. This qualitative observation is in accord with our earlier finding in the three flavour case (Eq. 10) that or slightly negative: in lattice units, the first order region tends to expand slightly as an imaginary chemical potential is turned on (see Fig. 7 (right)).

Similarly to the case, one expects the data along the whole line to be described by the leading term in the Taylor expansion,

(19) |

where now depends on . Conversion to continuum units proceeds as for , by determining the equivalent of Eqs. (8,12) for . We do this for fixed physical strange quark mass and scanning in , which now plays the role of the variable quark mass, and find . For the variation of the pseudo-critical coupling, Eq. (7), we obtain the coefficients given in Table 4, leading to and in this case.

5. | 1838(3) | 0. | 572(9) | — | 1. | 75(13) | 1. | 5 |
---|

These are similar in magnitude to the three-flavour values. Note that means that remains constant as the chemical potential is turned on. The decrease of the pseudo-critical temperature with real , given by , is then the dominant effect. It dictates that decreases when a real chemical potential is turned on. In other words, the first order region shrinks.

While our data for and small quark masses have larger errors which do not yet allow to constrain quantitatively, Fig. 9 and Eq. (12) leave little doubt that this coefficient is going to be negative along the whole upper part of the line. We thus arrive at the conclusion that, for the lattice spacing considered here, the curvature of the critical surface at is negative, and the first order region is shrinking when a real chemical potential switched on.

### 5.4 An alternative scenario for the QCD phase diagram

Let us take our results at face value for a moment and consider the implications if such a qualitative result also holds in the continuum limit. This leads to a scenario for the phase diagram which is at odds with common expectations. We find that the first order region in a plane of constant is actually shrinking with growing real . If the physical point is in the crossover region at , then switching on a chemical potential will not lead to an intersection with the critical surface as long as the latter is well described by its curvature at , i.e. for , or equivalently MeV. In the absence of any additional (and so far unknown) critical structure, there would then be no critical point or first order phase transition at all. The phase diagram of physical QCD would then only have the possible transition line separating the superconducting phase from nuclear matter, as in Fig. 11 (right).

Note that this scenario is perfectly consistent with all universality arguments and the known results for . This can be illustrated in the three flavour theory by considering the change of the -diagram with quark mass, as depicted in Fig. 12. All boundary conditions are met, in particular there is a first order phase transition at for quark masses smaller than . However, according to the negative sign for in Eqs. (17,9), the critical endpoint is now moving to the left with growing quark masses, until it disappears entirely.

It is natural to ask how reliable this unexpected scenario is regarding systematic errors. We have seen in Sec. 5.2 that lattices are very coarse, and we have discussed the enormous sensitivity of the critical values of the mass parameters to cut-off effects [32]. It is therefore expected that the values for shift once this analysis is repeated on finer lattices and/or with improved actions. The crucial question is whether the continuum extrapolated slope in Fig. 7(right) will be positive or negative and how it will balance out against the contribution from the pseudo-critical temperature. Or, equivalently, how the ordering of the zero and finite critical lines in Fig. 9 turns out in the continuum limit once physical units are used. In particular, our continuum conversion is sensitive to the non-perturbative beta function. A look at Eq. (12) gives a quick estimate of what is needed for a positive . Either would have to be positive of the order , or has to grow by a factor larger than 10 on the way to the continuum. Thus, while we presently cannot rule out that the picture reverts back to the standard scenario in the continuum limit, the opposite is obviously also possible as suggested by our data.

Finally, all of our arguments are based on the simplest scenario, in which the finite critical point of physical QCD is continuously connected to a critical point for some other mass values at . We cannot exclude a more complicated possibility, where the phase boundary of the color superconducting phase (see Fig. 12) would distort as a function of the quark masses, and give birth to a critical point distinct from the one we study, and which would survive for physical quark masses. This would correspond to a scenario with an additional critical surface in Fig.11 (left) above the one we studied here. We thus conclude that even the qualitative features of the QCD phase diagram cannot be regarded as settled yet.

### 5.5 The finely tuned critical end point

While based on our present data we are unable to make a reliable quantitative prediction for the location of the critical point for physical QCD, we have obtained important qualitative information regarding its quark mass dependence. Irrespective of the continuum extrapolated sign, all our evidence is that the absolute value of in Eq. (2) will be as naturally expected, while the effect of subleading terms is small up to MeV. This means that the critical quark masses are very weakly varying functions of the chemical potential, which is in line with the corresponding behaviour of the pseudo-critical temperature or, indeed, the equation of state. Consequently the inverse function is very strongly varying with quark mass, i.e. the finite critical surface emerges very steeply from the critical line in Fig. 2. Hence, even if in the continuum limit the conventional scenario with positive curvature is realized, the precise location of the critical end point will be exceedingly quark mass sensitive. A simple estimate using shows that, in order to have MeV, the physical point has to lie within of the chiral critical line, i.e. the physical quark masses would be fine tuned. While there is nothing forbidding such a situation, it appears rather unnatural. Moreover, if realized in nature, it would make a quantitative determination of through simulations exceedingly difficult. (For example, one might even need to treat the - and -quarks as non-degenerate).

## 6 Discussion and conclusions

As we have demonstrated in the preceding sections, step size errors are eliminated from our calculations. Finite volume effects are under control in the theory and for in , where we have performed our finite volume check. Note that this corresponds to pion masses larger than physical. For the point on the critical line with the lightest quark mass considered, we only have , and finite volume effects are to be expected. Ideally the part of the critical line in the neighbourhood of the physical point should be checked on a larger lattice.

However, by far the largest source of uncertainty is due to the coarse lattice spacing fm, as evidenced by several aspects of these computations. Strong cut-off effects reveal themselves when attempting to set a physical scale for the problem. Moreover, a change in the discretization of the Dirac operator on a lattice this coarse can change the pion mass corresponding to the second-order transition by a factor [7, 32] at . Finally, it has recently been pointed out that staggered simulations at finite suffer from additional discretization errors compared to [33] when , due to the eigenvalue structure when taking the fourth root of the determinant. For simulations at imaginary , the eigenvalues are pure imaginary, and this additional error is of , with possibly a large coefficient. A safe strategy thus is to first take the continuum limit of imaginary simulations, and then continue to real . For reweighting approaches at real one even expects errors.

In interpreting our findings and comparing with other work, it is important to take systematic uncertainties into account.
Given the cut-off effects, the sensitivity of the critical point to step size errors and, most notably, to the quark mass, it is clear that the discrepancy between our findings and those of [14] is nothing remarkably unusual, but merely reflects the large and different systematic uncertainties afflicting
these calculations.
In particular, in [14] the quark masses were held fixed in lattice units while was
increased.
Equivalently, was kept fixed. However, decreases under the
influence of a chemical potential, in a manner similar to Eq. 8,
so that the quark masses at the critical endpoint in [14] are about 5-10% smaller than physical.
This small deviation from a line of constant physics has a large impact on the
location of the critical point, because of the high sensitivity of the latter on
quark masses^{3}^{3}3In our approach, the conversion from lattice units
Eq. (10) to physical units Eq. (9) changes the coefficient ,
which is nearly zero, to , which is negative..
The effect is to artificially move the critical point to smaller
chemical potentials. The shifted masses may even reside in the first order
region, causing a critical point to be found even if in fact there is none,
consistently with the scenario discussed here.

Our study of the chiral critical surface also suggests that one cannot draw conclusions for physical QCD from
simulations of the critical point in QCD [26, 21].
This becomes clear when considering Fig. 1,
which describes the expectations.
If one moves from the physical point upwards by increasing to infinity,
the distance to the critical line increases considerably.
Given the high sensitivity of the critical chemical potential to this
distance (i.e. the small curvature of the chiral critical surface in
Fig. 2) which we observe, one should expect large differences
for between the and theories
^{4}^{4}4Moreover, for one expects a tricritical point
, with mean-field critical exponents which govern
the analytic form of the critical line , in contrast
to the case..

Resolving these various systematic issues and deciding which scenario for the phase diagram is realized in nature thus urgently requires further investigations of the theory with exact algorithms on finer lattices. Among the various finite approaches, our imaginary simulations require comparatively moderate computer resources to achieve this goal.

Acknowledgements: We are grateful to K. Rummukainen, Y. Shamir and B. Sve-titsky for discussions. We thank M. Stephanov for a thorough reading and many discussions. We also thank the Minnesota Supercomputer Institute, the RCNP Osaka and the HLRS Stuttgart for computer resources.

## References

- [1] E. Laermann and O. Philipsen, Ann. Rev. Nucl. Part. Sci. 53 (2003) 163 [arXiv:hep-ph/0303042].
- [2] R. D. Pisarski and F. Wilczek, Phys. Rev. D 29 (1984) 338; F. Wilczek, Int. J. Mod. Phys. A 7(1992) 3911 [Erratum-ibid. A 7 (1992) 6951]; K. Rajagopal and F. Wilczek, Nucl. Phys. B 399 (1993) 395 [arXiv:hep-ph/9210253].
- [3] R. V. Gavai, J. Potvin and S. Sanielevici, Phys. Rev. Lett. 58 (1987) 2519; J. B. Kogut and D. K. Sinclair, Nucl. Phys. B 295, 480 (1988); F. R. Brown et al., Phys. Rev. Lett. 65, 2491 (1990).
- [4] S. A. Gottlieb, J. Kuti, D. Toussaint, A. D. Kennedy, S. Meyer, B. J. Pendleton and R. L. Sugar, Phys. Rev. Lett. 55 (1985) 1958.
- [5] P. Hasenfratz, F. Karsch and I. O. Stamatescu, Phys. Lett. B 133 (1983) 221.
- [6] C. Alexandrou, A. Borici, A. Feo, P. de Forcrand, A. Galli, F. Jegerlehner and T. Takaishi, Phys. Rev. D 60 (1999) 034504 [arXiv:hep-lat/9811028]; A. Dumitru, D. Roder and J. Ruppert, Phys. Rev. D 70 (2004) 074001 [arXiv:hep-ph/0311119].
- [7] F. Karsch, E. Laermann and C. Schmidt, Phys. Lett. B 520 (2001) 41 [arXiv:hep-lat/0107020].
- [8] N. H. Christ and X. Liao, Nucl. Phys. Proc. Suppl. 119 (2003) 514.
- [9] P. de Forcrand and O. Philipsen, Nucl. Phys. B 673 (2003) 170 [arXiv:hep-lat/0307020].
- [10] S. A. Gottlieb, W. Liu, D. Toussaint, R. L. Renken and R. L. Sugar, Phys. Rev. D 35 (1987) 2531.
- [11] M. A. Clark, A. D. Kennedy and Z. Sroczynski, Nucl. Phys. Proc. Suppl. 140 (2005) 835 [arXiv:hep-lat/0409133].
- [12] T. Herpay and Z. Szep, arXiv:hep-ph/0604086; T. Herpay, A. Patkos, Z. Szep and P. Szepfalusy, Phys. Rev. D 71 (2005) 125017 [arXiv:hep-ph/0504167]; J. T. Lenaghan, Phys. Rev. D 63 (2001) 037901 [arXiv:hep-ph/0005330]; Y. Hatta and T. Ikeda, Phys. Rev. D 67 (2003) 014028 [arXiv:hep-ph/0210284].
- [13] Z. Szep, PoS JHW2005 (2006) 017 [arXiv:hep-ph/0512241].
- [14] Z. Fodor and S. D. Katz, JHEP 0404 (2004) 050 [arXiv:hep-lat/0402006].
- [15] O. Philipsen, PoS LAT2005 (2005) 016 [PoS JHW2005 (2006) 012] [arXiv:hep-lat/0510077].
- [16] Ph. de Forcrand and O. Philipsen, Nucl. Phys. B 642 (2002) 290 [arXiv:hep-lat/0205016].
- [17] A. Roberge and N. Weiss, Nucl. Phys. B 275 (1986) 734.
- [18] M. D’Elia and M. P. Lombardo, Phys. Rev. D 67, 014505 (2003) [arXiv:hep-lat/0209146].
- [19] A. Hart, M. Laine and O. Philipsen, Phys. Lett. B 505 (2001) 141 [arXiv:hep-lat/0010008].
- [20] S. Choe et al., Phys. Rev. D 65 (2002) 054501.
- [21] C. R. Allton et al., Phys. Rev. D 71 (2005) 054508 [arXiv:hep-lat/0501030]; C. R. Allton, S. Ejiri, S. J. Hands, O. Kaczmarek, F. Karsch, E. Laermann and C. Schmidt, Phys. Rev. D 68 (2003) 014507 [arXiv:hep-lat/0305007].
- [22] K. Binder, Z. Phys. B 43 (1981) 119.
- [23] J. B. Kogut and D. K. Sinclair, arXiv:hep-lat/0504003.
- [24] M. A. Clark, P. de Forcrand and A. D. Kennedy, PoS LAT2005 (2005) 115 [arXiv:hep-lat/0510004].
- [25] A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. 63 (1989) 1195.
- [26] R. V. Gavai and S. Gupta, Phys. Rev. D 71 (2005) 114014 [arXiv:hep-lat/0412035].
- [27] S. A. Gottlieb, Nucl. Phys. Proc. Suppl. 20 (1991) 247; C. R. Allton et al., Phys. Rev. D 66, 074507 (2002) [arXiv:hep-lat/0204010].
- [28] D. K. Sinclair and J. B. Kogut, arXiv:hep-lat/0609041.
- [29] Ph. de Forcrand, S. Kim, S. Kratochvila, T. Takaishi, arXiv:hep-lat/0510069.
- [30] M. D’Elia, A. Di Giacomo and C. Pica, Phys. Rev. D 72 (2005) 114510 [arXiv:hep-lat/0503030].
- [31] J. B. Kogut and D. K. Sinclair, Phys. Rev. D 73 (2006) 074512 [arXiv:hep-lat/0603021].
- [32] F. Karsch, C. R. Allton, S. Ejiri, S. J. Hands, O. Kaczmarek, E. Laermann and C. Schmidt, Nucl. Phys. Proc. Suppl. 129 (2004) 614 [arXiv:hep-lat/0309116].
- [33] M. Golterman, Y. Shamir and B. Svetitsky, arXiv:hep-lat/0602026.