A finite-density transition line for QCD with 695 MeV dynamical fermions

We apply the relative weights method to SU(3) gauge theory with staggered fermions of mass 695 MeV at a set of temperatures in the range $151 \le T \le 267$ MeV, to obtain an effective Polyakov line action at each temperature. We then apply a mean-field method to search for phase transitions in the effective theory at finite densities. The result is a transition line in the plane of temperature and chemical potential, with an endpoint at high temperature, as expected, but also a second endpoint at a lower temperature. We cannot rule out the possibilities that a transition line reappears at temperatures lower than the range investigated, or that the second endpoint is absent for light quarks.


I. INTRODUCTION
The effective Polyakov line action (PLA) of a lattice gauge theory is the theory which results from integrating out all of the degrees of freedom of the theory, subject to the condition that the Polyakov lines are held fixed, and it is hoped that this effective theory is more tractable than the underlying lattice gauge theory (LGT) when confronting the sign problem at finite density. The general idea was pioneered in [1], and the derivation of the PLA from the underlying LGT has been pursued by various methods, e.g. [2][3][4]. The relative weights method [5] is a simple numerical technique for finding the derivative of the PLA in any direction in the space of Polyakov line holonomies. 1 Given some ansatz for the PLA, depending on some set of parameters, we can use the relative weights method to determine those parameters. Then, given the PLA at some fixed temperature T , we can apply a mean field method to search for phase transitions at finite chemical potential µ. This is the strategy which we have outlined in some detail in [6], where some preliminary results for finite densities were presented. The relative weights method has strengths and weaknesses; on the positive side the approach is not tied to either a strong coupling or hopping parameter expansion, and the non-holomorphic character of the fermion action is irrelevant. The main weakness is that the validity of the results depends on a good choice of ansatz for the PLA. We have suggested, for exploratory work, an ansatz for the PLA inspired first by the success of the relative weights method applied to pure gauge theories [5], and secondly by the form of the PLA obtained for heavy-dense quarks.
In this article we follow up on the work in ref. [6] to obtain a tentative transition line in the µ − T plane for SU (3) gauge theory with dynamical staggered unrooted quarks of mass 695 MeV. It is generally believed that this line has an endpoint at high temperatures, and this is what we find. A second, unexpected finding is that there is also an endpoint at a lower temperature. 2 Whether a transition line reappears at still lower temperatures, outside the range we have investigated, or whether the second transition point disappears for lighter quarks, or whether this second transition point is instead indicative of some deficiency in our ansatz for the PLA, remains to be seen.
In the next section we briefly review the relative weights method and associated mean field technique at finite density, referring to our previous work [6] for some of the technical details. Section 3 contains our results, followed by conclusions in section 4. A result for a recently introduced observable ξ /ξ 2nd [10], which is sensitive to excitations above the lowest lying excitation, is presented in an appendix.

II. RELATIVE WEIGHTS
It is simplest to work in temporal gauge, where we can fix the timelike links to the identity everywhere except on one timeslice, say at t = 0, on the periodic lattice. In this gauge the timelike links at t = 0 are the Polyakov line holonomies, which are held fixed, and the PLA S P for an SU(N) gauge theory is defined by where S L is the lattice action for an SU(N) gauge theory with dynamical fermions, and we also define the Polyakov line P x x x = 1 N TrU x x x . In this article we use the standard SU(3) Wilson action for the gauge field, and unrooted staggered fermions as the dynamical matter fields. Given the PLA at µ = 0, the PLA at finite µ is obtained by the simple replacement where N t is the lattice extension in the time direction, with µ in lattice units. Let us consider some path through the space of Polyakov line holonomies U x x x (λ ) parametrized by λ . The relative weights method allows us to compute the derivative dS P /dλ anywhere along the path. Let U ′ x x x and U ′′ x x x represent the Polyakov line holonomy field at parameter λ 0 ± 1 2 ∆λ respectively. Defining S ′ L , S ′′ L as the lattice actions with the Polyakov line holonomies fixed to U ′ x x x ,U ′′ x x x respectively, and , then in temporal gauge we have by definition where ... ′′ indicates that the expectation value is to be taken in the probability measure The expectation value in the last line of (3) can be calculated by standard lattice Monte Carlo, only holding fixed timelike links at t = 0. From this calculation we find the derivative dS P /dλ ≈ ∆S P /∆λ . The SU(2) and SU(3) gauge groups are special in the sense that the trace P x x x determines the holonomy U x x x . So we expand and compute the derivatives ∂ S P /∂ a k k k with respect to a set of specific Fourier components a k k k , keeping the remaining components fixed to values taken from a thermalized configuration. Our ansatz for the PLA is where p = 1 for unrooted staggered fermions, and p = 2N f for Wilson fermions, where N f is the number of flavors. The determinant factors are motivated by the PLA for heavy dense quarks [1,11] in which h = (2κ) N t , with κ the hopping parameter for Wilson fermions, or κ = 1/2m for staggered fermions. In our ansatz, h becomes a fit parameter. The part of the action involving the kernel K(x x x−y y y) is motivated by previous successful treatments [5] of pure gauge theory and gauge-Higgs theory. All in all, using where K(x x x − y y y) = 1 on an L 3 three-space volume, and choosing p = 1, we have The Fourier transform K(k k k) of K(x x x − y y y) is determined numerically at µ = 0. For k k k = 0, where If h ≪ 1, then dropping terms of O(h 2 ) and higher the derivative simplifies to where the "R" superscript refers to the real part. For k k k = 0, again dropping terms of order h 2 and higher, The derivatives on left-hand sides of (14) and (15) are determined by the method of relative weights at a variety of α. By plotting those results vs. α and fitting the data to a straight line, K(k k k) is determined from the slope. The h parameter can in principle be determined by the y-intercept of the K(0) vs. α data. A better method, for reasons to be discussed below, is to choose h by requiring that P computed in the PLA at µ = 0 agrees with P computed in the underlying LGT,

III. DERIVING THE POLYAKOV LOOP ACTION
It is clear, since there should be only finite range correlations in the PLA, that K(R) must die off faster than a power at large R, since otherwise one could bring down terms in the effective action to produce a power-law falloff in the Polyakov line correlator. However, since K(k k k) is determined at only a small set of k k k values on a 16 3 volume, it is necessary to fit the data to some analytic expression in order to carry out the inverse Fourier transform to K(R). It is of interest to see whether there are indications of the required rapid falloff at large R. As in our previous work [5,6], we fit the data for K(k k k) at k k k = 0 by a double straight-line fit where is the lattice momentum, and carry out the inverse Fourier transform, but taking K(0) from the data rather than the fit. The double-line fit at β = 5.7, ma = 0.6, N t = 6 is shown in Fig. 1, and the corresponding inverse Fourier transform to K(R) is displayed in Fig. 2. Qualitatively, the behavior shown is typical of the data for K(R) in all of our simulations. What we observe is that K(R) falls off roughly as 1/R 4 up to R ≈ 4, as seen from a straight-line fit on the log-log plot of Fig. 2(a). Beyond R = 4 K(R) passes through zero ( Fig. 2(b)), after which the data oscillates and is rather noisy. We believe that the noisy behavior is an artifact of having a small lattice, and a relatively small number of data points for K(k k k). But we do see, up to the onset of noise, precisely the rapid drop to zero which is expected on general grounds. We therefore discard the values of K(R) beyond R = R cut , and reset those values to where R cut is taken as the point where K(R) is either zero, or else reaches a minimum near zero before oscillating.
We have also tried to fit the data for K(k) to a third-order polynomial. Although the resulting K(R) agrees with our double-line method up to R ≈ 4, it then goes asymptotically to a non-zero constant, which is not the right behavior. We therefore prefer the original double line fit used in our previ- ous work, which at least indicates that K(R) should be negligible beyond some R cut , and which also gives us a fairly precise criterion for the choice of R cut . In fact by this criterion we find R cut = 4.6 in all cases. We report on the effects of small variations in the choice of R cut below in section V.
With K(R) determined by the procedure just outlined, we find h by simulating the PLA in eq. (11) at µ = 0 with a series of trial values of the h parameter, until P computed in the PLA agrees with the value computed in the LGT up to error bars. This is not entirely straightforward. As we found in our earlier work [6], the highly non-local term containing K(x x x − y y y) in the action leads to metastable states which, in a lattice simulation of the PLA, depend on the initialization and persist for thousands of iterations. A cold start, which sets P(x x x) = 1 initially, generally leads to the system in the "deconfined" phase, with a large expectation value P , even at h = 0. This is in strong disagreement with the LGT. Instead we initialize the system at P(x x x) = 0, which then stays in the "confined" phase.
In principle h could be determined from our data via eq. (12), or the simplified version in eq. (14). We have not followed this approach for two reasons. First, the value of h turns out to be very small at β = 5.7 and below, and the value determined from a straight-line fit of data for the left hand side of (14) vs. α has a very large error bar. As an example, we plot this data, and the corresponding straight-line fit, in Fig.  3(a) for the case of β = 5.7. In this case a straight-line fit gives h = 0.0056 ± 0.0023, where the error bar is about 40% of the estimate. The alternative procedure used in this paper, of choosing h to get agreement for the Polyakov line value in the PLA and the LGT at µ = 0, gives h = 0.0042, which is actually well within the sizeable error bar associated with the straight-line fit procedure. But the relative error ∆h/h only gets worse at smaller β , and for this reason we have resorted to the Polyakov-fit method to determine h. h rises rapidly above β = 5.7, but here we encounter a different kind of problem: the data for dS/da 0 as a function of α does not fit a straight line, particularly for α at the Polyakov line expectation value P . This is illustrated in Fig. 3(b). One could try to fit the right hand side of (14) to the tangent at α = P , but we believe it is better to use a consistent procedure for all β values, so we have derived h above and below β = 5.7 in the same way, choosing an h to get P in agreement with the lattice gauge theory value.
Having obtained h, we can then go on to compute the Polyakov line correlator at µ = 0 The PLA still has a sign problem at µ = 0, which we deal with via a mean field method [12]. We summarize here only the essential points; a detailed derivation may be found in the cited references.
Starting from the partition function of the effective theory where D x x x (µ, TrU, TrU † ) is the product of the determinant factors (7) and (8), we rewrite where V = L 3 is the lattice volume, and we have defined Parameters u and v are to be chosen such that E 0 can be treated as a perturbation, to be ignored as a first approximation. In particular, E 0 = 0 when and these conditions turn out to be equivalent to the stationarity of the mean field free energy. The mean field approximation is obtained, at leading order, by dropping E 0 , in which case the partition function factorizes, and can be solved analytically as a function of u and v. After some manipulations (cf. [12]), one finds the mean field approximations u, v to TrU x and TrU † x respectively, by solving the pair of equations where A = J 0 v, B = J 0 u. The expression G(A, B) is given by where D −s i j is the i, j-th component of a matrix of differential operators The mean field free energy density f m f and fermion number density n are The stationarity conditions (24) may have more than one solution, and here it is important to take account of the existence of very long lived metastable states in the PLA. The state at µ = 0 which corresponds to the LGT is the one obtained by initializing at P x x x = 0, and in the mean field analysis this is actually not the state of lowest free energy (its stability in a Monte Carlo simulation is no doubt related to the highly nonlocal couplings in the PLA). By analogy, at finite µ we look for solutions of (24) by starting the search at u = v = 0, regardless of whether another solution exists at a slightly lower f m f . The mean field calculation at any chemical potential µ requires three numbers: a 0 and J 0 , which are both derived from K(x x x), and h. Denote by h(mfd) the value of h for which mean field gives a Polyakov line P , at µ = 0, in agreement with lattice gauge theory. This can be compared to the value of h, denoted h(PLA), for which the Polyakov line computed in the PLA at µ = 0 agrees with the lattice gauge theory value. The values of h(mfd) and h(PLA) at the points we have computed can be compared in Table I; in general they are not far off, which is some evidence of the validity of the mean field approach in this application.

Our
LGT numerical simulations were carried out in SU (3) We take the quark mass in lattice units to be ma = 0.6 at β = 5.7. This corresponds to a mass of m = 695 MeV, and temperature T = 193 MeV in physical units. We keep the physical mass and the extension N t = 6 in the time direction fixed, and vary the temperature by varying the lattice spacing, i.e. by varying β . Given the PLA at µ = 0, the Polyakov line expectation values TrU x x x and TrU † x x x at finite µ are calculated by the mean field method outlined above, with a sample of our results displayed in Fig. 5. A discontinuity in a plot of TrU x x x , TrU † x x x vs. µ is the sign of a transition at finite density, and conversely the absence of any discontinuity indicates the absence of any transition. When a transition occurs at some value of chemical potential µ 1 , then there is a second transition at some µ 2 > µ 1 . However, while the first transition occurs at some relatively low density (in lattice units) on the order of n ≈ 0.2, the second transition always occurs at a density close to the saturation value, which for staggered fermions is n = 3. An example of density n vs. µ at β = 5.75 is shown in Fig. 6; transitions occur at the sharp jumps in density. Since the saturation value is a lattice artifact, we do not attach much physical significance to the second transition at µ = µ 2 . Transition points listed in our tables and plots all correspond to µ = µ 1 .
Our full set of results, including the parameters used in the mean field calculation, is shown in Table I. As a check on the mean field calculation we have calculated h in two ways: • by choosing h such that the mean field calculation at µ = 0 reproduces the Polyakov line value of lattice gauge theory. These values are denoted h(mfd) in Table  I. The corresponding first transition µ 1 /T , and µ 1 in MeV, are displayed in the columns adjacent to h(mfd).
• by choosing h such that a simulation of the PLA at µ = 0 reproduces the Polyakov line value of lattice gauge theory, as described in section III. These values are denoted h(PLA) in Table I. The corresponding first transition µ 1 /T , and µ 1 in MeV, are displayed in the columns adjacent to h(PLA).
The fairly close agreement between h(mfd) and h(PLA), in what amounts to a spin system in D = 3 dimensions, can be attributed to the highly non-local nature of the PLA. Mean field treatments work best when each degree of freedom is coupled to many other degrees of freedom. For nearest-neighbor couplings, this generally means that the treatment is only accurate in higher dimensions. But it seems that the relatively long range of K(R), which results in each Polyakov line being coupled to many other lines, serves the same purpose as high dimensions in a nearest-neighbor theory. We should also note the good agreement found in ref. [12] between Langevin simulations at finite density in a spin system of this type, with the mean field treatment we have discussed.
Our results for the µ, T transition points (using h=h(PLA)) for staggered, unrooted quarks of mass 695 MeV, are plotted in Fig. 7(a). This figure is the main result of our paper, and holds for the temperature range 151 ≤ T ≤ 267 MeV that we have investigated. We see that the phase transition line exists to an upper temperature of T ≈ 241 MeV, where there is a critical endpoint. The fact that there is a critical endpoint at high temperature was expected. What was unexpected is that the transition points seem to disappear at lower temperatures. The solid line is a comparison to the analytic continuation expression of d'Elia and Lombardo [14], to be discussed further below.
There are two transition points, computed at β = 5.65 and β = 5.66, which appear to be outliers on this plot, in the sense that they lie quite far from the analytic continuation curve. These β values are in fact the smallest couplings at which we still see a transition, and differ sharply from the transition point at the next lowest β value, at β = 5.68 with the transition at (µ, T ) = (961, 184) MeV. We suspect that these two points are indicative of some instability in the mean field calculation in the immediate neighborhood of the β value where the phase transitions disappear, and we are inclined to discount them. Taken literally, they would suggest that the transition line moves to lower values of µ as T decreases below 184 MeV, which seems unlikely. One indication that something may be going wrong at these lowest β values is the fact that at β = 5.66, unlike at all our other data points, there is a very substantial disagreement between the transition point at µ = 928, which is derived using h(mfd), and the transition point µ = 849 MeV, derived using h(PLA). A comparison of transition points derived in the mean field approach  Note that the transition at higher chemical potential occurs close to the saturation value, which is n = 3 for staggered fermions. Since saturation is a lattice artifact, this second transition is probably unphysical.
using h(mfd) and h(PLA) is displayed in Fig. 8. In this figure we have drawn a short horizontal line connecting transition points at β = 5.66, to indicate their separation. A possible source of systematic error could lie in our procedure for choosing R cut , and it is worthwhile to check the sensitivity of our results to a small variation in that parameter. This is shown in Table II, where we have varied R cut (in lattice units) by ±0.2. These should not be regarded as error bars, exactly, since the variation ±0.2 is rather arbitrary, but only as some indication of the sensitivity of our results to the value of R cut .
In any case, our results do indicate that the transition line does not continue all the way to T = 0, and seems to end at about T = 184, µ = 961 MeV. We cannot rule out the possibility that a high density transition reappears at some temperature lower than the lowest temperature (151 MeV) that we have considered. Or perhaps the second critical endpoint goes away for light quarks. These possibilities we reserve for later investigation.

A. Comparison with other results
There do exist results for the critical line in the µ − T plane, expected to be valid at small µ, that have been obtained by other methods. Because these usually involve a choice of   lattice fermions (Wilson, staggered), number of flavors, and quark masses different from our own, it is a little difficult to make a direct comparison with our data. Nevertheless, it is interesting make such comparisons anyway, for whatever they may be worth.
Perhaps the most relevant comparison is to the analytic continuation method of d'Elia and Lombardo [14], who, like us, work with four flavors of staggered fermions. Their approach is to find the transition line for imaginary chemical potential, and fit to a polynomial which is then analytically continued to real chemical potential. The result (like the related Taylor expansion approach), is only believable at small µ. However, these authors set ma = 0.05, which means that their physical mass changes at different β . What is found is that where T c is the critical (or crossover) temperature at µ = 0. In order to make a comparison with our work we just take T c as a fitting parameter, with the expectation that it can't be too far from the quenched case of around T c = 250 MeV (in fact the fit comes out to be T c ≈ 233 MeV). The comparison of (29) with our result is shown in Fig. 7(a), and there appears to be remarkably good agreement between the analytic continuation result, and the transition points we have found via relative weights, apart from the two outliers mentioned above.
We have also made the comparison to the results for heavy dense quarks obtained by Aarts et al. [15] via the complex Langevin method. That work employs Wilson fermions with two flavors, and a hopping parameter of κ = 0.04; i.e. extremely massive quarks. These authors fit their data for the critical temperature at each µ to a 2nd order polynomial and µ 0 = − ln(2κ). This choice of µ 0 is motivated by the hopping parameter expansion, and the "Silver Blaze" phenomenon (nothing happens at T = 0 until µ = µ 0 ). Aarts et al. [15] report the constants on their largest (10 3 ) lattice volume. There is some volume dependence in these constants, and the result is presumably only valid for small hopping parameter and large chemical potential. In order to make some kind of comparison with this with heavy-dense approach, we have just taken µ 0 to give the closest fit to our data points. The result of this fit is shown in Fig. 7(b). It's hard to know how seriously to take this comparison, since Aarts et al. are using two flavors of Wilson fermions, as opposed to our four staggered flavors, and are working with extremely heavy quarks. In any case, even allowing for a best fit value of µ 0 , the heavy dense result seems quite far from our data, and perhaps this is not very surprising.

VI. CONCLUSIONS
We have found a first-order phase transition line for SU(3) gauge theory with dynamical unrooted staggered fermions of mass 695 MeV, by the method of relative weights combined with mean field theory, in the plane of chemical potential µ and temperature T . The critical line lies in a finite temperature range between (µ, T ) = (175, 241) and (µ, T ) = (961, 184) MeV and the transition points lie, for the most part, along a line determined from analytic continuation of results at imaginary µ [14].
We offer this result with reservations. There is certainly a degree of arbitrariness in our approach, particularly in the choice of ansatz for the PLA. We have used a product of local determinants for the µ-dependent part of the PLA, and a non-local bilinear form for the µ-independent part. This form has given us an excellent match to the Polyakov line correlator at µ = 0, computed in the LGT. But there is no guarantee that a more complicated form, involving e.g. Polyakov lines in higher representations, trilinear couplings, etc., is not required at high densities. It would be interesting to probe the existence and possible importance of such terms, supplementing the method of relative weights and mean field with, e.g., the method of inverse Monte Carlo [4,16] and/or the approach of Bergner et al. [3]. It would also be interesting to see whether the expected transition line reappears, for m = 695 MeV, at temperatures below 151 MeV, or how the situation changes with lighter quark masses. We reserve these questions for later study. Gauge theories have a complicated spectrum, and correlators at intermediate distances should be sensitive to more than just the lowest-lying excitation. Caselle and Nada [10] have introduced an interesting observable which tests for the presence, in an effective Polyakov line action, of a spectrum of excitations contributing to two-point correlators. We briefly review their idea in this section, as implemented at µ = 0.
Denoting P x x x = P x,y,z , define the average of Polyakov lines in a plane on an L 3 lattice volume Let G ∞ (τ) be G(τ) in the L → ∞ limit. Following [10] we consider If we approximate the sums over τ by integrals, and assuming that G ∞ (τ) has the form (A.2) for all τ, we find that Now suppose that the sum is dominated by the first term, with all other terms negligible. In that case ξ defined in (A.3) On the lattice we have lattice periodicity, and we cannot take the L → ∞ limit to get G ∞ (τ). The proposal is instead to approximate (A.5) by where ξ is determined from (A.3).
There are, of course, possible sources of systematic error in the choice of τ max and the extraction of ξ . But for a PLA derived for a system of dynamical fermions there is an additional ambiguity. Polyakov lines at a separation less than the stringbreaking scale represent a quark-antiquark system joined by a flux tube, and the contributions to G(τ) come from a spectrum of string-like excitations. Beyond the string-breaking scale, the Polyakov lines represent two bound states, namely a static quark + light antiquark, and a static antiquark + light quark. In this case the contributions to the connected correlator are associated with hadron exchange. So in this case the ratio ξ /ξ 2nd is picking up contributions from quite different regimes. In view of this, we have calculated the ratio ξ /ξ 2nd at β = 5.65, where the expectation value of the Polyakov line is very small, and the contributions to ξ /ξ 2nd are coming mainly from the string-like spectrum, at least on the comparatively small lattice volume of 16 4 used in our simulations. We obtain ξ = 2.196(11) from the fit (A.3), with data and fit shown in Fig. 9 shown, and using (A.7) with τ max = 5 we finally ob-tain ξ ξ 2nd = 1.27 ± 0.03 (A. 8) which is comparable to some of the results for pure SU(2) lattice gauge theory quoted in ref. [10]. The fit to a single cor-relation length (A.3) is expected to be very accurate at large τ, and the influence of higher excitations would only be evident at smaller τ. Note that in Fig. 9 the fit (A.3) is in fact quite accurate for all τ > 0, so in this case the deviation of the ratio from unity must be attributed mainly to the data point at τ = 0. We must stress again that this result is obtained on a relatively small lattice volume, and is subject to the caveats already mentioned.