Preliminary QCD phase transition line for 695 MeV dynamical staggered fermions from effective Polyakov line actions

We present a phase diagram for SU(3) lattice gauge theory with 695 MeV dynamical staggered fermions, in the plane of temperature and chemical potential, obtained from effective Polyakov line actions. The derivation is via the method of relative weights, and the effective theories are solved at finite chemical potential by mean field theory.


Introduction
One of the most active areas in strong interaction physics concerns the behavior of QCD in extreme conditions, i.e., high temperature and/or high baryon density. At high temperatures we enter the realm of the quark gluon plasma, whose properties have been or will be probed by experiments at RHIC, the LHC, and the FAIR facility (now under construction). Not much is known for sure about hadronic matter at high baryon density. QCD is believed to have a complex phase structure in the temperature-density plane, and a number of exotic phases (quarkyonic, glasma, color-flavor locked superconductor...) have been suggested . One would especially like to know the position of the critical end point of the confinement-deconfinement transition. Many talks on the subject of QCD in extreme environments begin with a sketch of the possible phase diagram, but such sketches are, so far, all conjecture. Nobody knows whether these exotic phases really exist, or exactly where they are located in the temperature/density plane. So the first order of business, for a theorist, is to nail down the phase diagram. By far the most important tool in the investigation of non-perturbative QCD is the method of importance sampling in lattice Monte Carlo simulations. But when one attempts to apply this tool to study QCD at high baryon density, a serious obstacle − the "sign problem" − is encountered. Different strategies were explored, and the reviews at the yearly lattice conferences [1][2][3][4][5][6][7][8][9] summarize the progress. Finite densities are introduced in statistical systems via the introduction of a chemical potential, but when this standard method is applied in QCD, the fermion determinant becomes complex and cannot be interpreted as a probability measure. Then standard importance sampling, e.g., via the hybrid Monte Carlo method, breaks down completely, and some other method or methods must be devised to handle the problem of complex actions. The sign problem crops up not only in QCD at high baryon density, but also in various condensed matter systems involving highly correlated fermions. So apart from the general conceptual challenge, it would be rewarding in many areas of physics if we could figure out how to apply numerical methods to systems with Speaker, e-mail: hroman@kph.tuwien.ac.at complex actions. There are many promising avenues, of course, in particular the complex Langevin equation [10], Lefshetz thimbles [11], and the histogram approach [12]. However, no method so far is entirely free from objections, e.g., in the case of the complex Langevin equation, the problem is that the fermionic part of the action is non-holomorphic, due to the fact that the logarithm of the fermion determinant has a branch cut. It has been shown by Mollgaard and Splittorff [13], that this can lead complex Langevin evolution to incorrect results. As for the Lefshetz thimble approach, the rigorous statement is that all thimbles must be summed over, but the proposed application to QCD is to simulate quantum fluctuations around a single thimble. It is not yet clear that this is the right thing to do, and in simple examples the restriction to a single thimble gives the wrong answer. Finally, the histogram approach was studied critically by Joyce Myers, Kim Splittorff, and Jeff Greensite [14,15]. What was found is that the fundamental assumption of the histogram approach, namely, that the phase of the fermion determinant has a Gaussian distribution, is probably not right, and that even tiny deviations from a Gaussian distribution, on the order of inverse powers of the lattice volume, are sufficient to invalidate conclusions based on a Gaussian distribution alone.
Our approach to the sign problem in QCD is to map QCD with a chemical potential into a simpler effective theory, namely, the effective Polyakov line action (henceforth "PLA"), via the relative weights method and then deal with the sign problem via mean field theory, which is a surprisingly accurate method for solving effective actions of this kind [14]. The general idea was pioneered in [16], and the derivation of the PLA from the underlying LGT has been pursued by various methods, e.g. [17][18][19][20][21][22][23]. The relative weights method [24,25] is a simple numerical technique for finding the derivative of the PLA in any direction in the space of Polyakov line holonomies. 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 µ. The phase diagram of the effective theory will mirror the phase diagram of the underlying gauge theory. The method was successfully tested in SU(2) and SU(3) pure gauge and gauge-Higgs theories [24][25][26] and first results with dynamical fermions were presented in [27][28][29][30]. Here we follow up on this work and 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 129 ≤ T ≤ 260, 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. We obtain a tentative transition line in the µ − T plane that has an endpoint at high temperatures, as expected, and a second, unexpected endpoint at a lower temperature. It remains to be seen, whether a transition line reappears at still lower temperatures, 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. See also [31] for more details.

Polyakov Line Action and Relative Weights
The effective Polyakov line action S P is the theory obtained by integrating out all degrees of freedom of the lattice gauge theory, under the constraint that the Polyakov line holonomies are held fixed. It is convenient to implement this constraint in temporal gauge (U 0 (x, t 0) = ), so that where φ denotes any matter fields, scalar or fermionic, coupled to the gauge field, and S L is the SU(3) lattice action (note that we adopt a sign convention for the Euclidean action such that the Boltzmann weight is proportional to exp[+S ]). To all orders in a strong-coupling/hopping parameter expansion, the relationship between the PLA at zero chemical potential µ = 0, and the PLA corresponding to a lattice gauge theory at finite chemical potential, is given by So the immediate problem is to determine the PLA at µ = 0. Let us define the Polyakov line in an SU(N) theory to refer to the trace of the Polyakov line holonomy P Relative weights enables us to compute the derivative of the effective action S P along any path with spatial lattice extension L = 16. The procedure is to run a standard Monte Carlo simulation, generate a configuration of Polyakov line holonomies U x , and compute the Polyakov lines P x = TrU x . We then set the momentum mode a k = 0 in this configuration to zero, resulting in the modified where f is a constant close to one ( f = 1 is only possible in the large volume, α → 0 limit). From the holonomy configurations U x , U x we compute derivatives of S P with respect to the real part a R k of the The expectation value is straightforward to compute numerically, by fixing the Polyakov holonomies and calculating the action differences, and from the logarithm we determine ∆S P . Our proposal is to fit the relative weights data to an ansatz for S P based on the massive quark effective action [16,[32][33][34][35] where both the kernel K( x − y) and the parameter h are to be determined by the relative weights method. The full action is surely more complicated than this ansatz; the assumption is that these terms in the action are dominant, and the effect of a lighter quark mass is mainly absorbed into the parameter h and kernel K( x − y). We then have the derivative of the action with respect to momentum modes a k of the Polyakov lines The left-hand side is computed via relative weights at a variety of α = 0.01, . . . , 0.06, and plotting those results vs. α, K( k) is determined from the slope, while h is given by the intercept at α = 0 from the zero mode k = 0, since the second term basically vanishes for k 0, see Figure 1. We fine-tune h by requiring that P agrees in the PLA with the underlying LGT. From the kernel K( k) in momentum space we derive the kernel K( x − y) via inverse Fourier transformation and find that it behaves like K(r) = c/r 4 , see Figure 2a. At finite lattice spacing it deviates from the simple ansatz for large r, hence we introduce a cutoff R cut and set the kernel to zero for r > R cut in our PLA simulations. A first consistency check of our PLA simulations is that we reproduce the correct Polyakov line correlator of the underlying lattice gauge theory, as shown in Figure 2b for β = 5.7.

Finite Density and Mean Field Theory
For µ 0 the effective PLA has still a sign problem, which we solve via mean field theory as discussed in [36,37]. The treatment of S U(3) spin models at finite µ is just a minor variation of standard a) mean field theory at zero chemical potential. The basic idea that each spin is effectively coupled to the average spin on the lattice, not just nearest neighbors, is favored by the effective PLA with the non-local kernel K( x − y). We introduce two magnetizations u and v for TrU x and TrU † x which are determined by minimizing the free energy. Therefore we rewrite the kernel part of the action (7) to separate local and non-local terms where V = L 3 is the lattice volume, and we have defined u and v are chosen such that E 0 can be treated as a perturbation, in particular, E 0 = 0 when u = TrU x , v = TrU † x , 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. [36,37]), 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 G(A, B) and The stationarity conditions (11) 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 = 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 non-local couplings in the PLA). By analogy, at finite µ we look for solutions of (11) by starting the search at u = v = 0, regardless of whether another solution exists at a slightly lower f m f . A first check of the mean field approach is that we reproduce the correct expectation value of the Polyakov loop from LGT at µ = 0.

The QCD Phase Transition Line
We derive effective Polyakov line actions from lattice simulations of SU(3) Wilson gauge action and dynamical staggered fermions on 16 3 × 6 lattices, for a variety of temperatures and lattice masses corresponding to a physical quark mass m q = 695MeV, as summarized in Table 1. To set the scale we use the lattice spacing a from the Necco-Sommer formula [38]. 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 β. Polyakov loop correlators in the LGT and the PLA at µ = 0 are in very good agreement in all cases as well as the Polyakov line expectation values after mean field treatment.
Turning on the chemical potential µ, TrU x and TrU † x are calculated by the mean field method outlined above, with a sample of our results displayed in Fig. 3. A discontinuity 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.1, the second transition always occurs at a density close to the saturation value, which for staggered fermions is n = 3. Since the saturation value is a lattice artifact, we do not attach much physical significance to the second transition.
From the µ 1 vs. T data shown in the table, we plot a transition line in the µ − T plane for staggered, unrooted quarks of mass 695 MeV, Fig. 4b. This figure is the main result of our paper, and holds for the temperature range 129 ≤ T ≤ 260 MeV that we have investigated. We see that the phase transition line exists to an upper temperature of T ≈ 220 MeV, where there is a critical endpoint. The fact that there is a critical endpoint at high temperature was expected. What was unexpected is the existence of a second critical endpoint at a lower temperature of T ≈ 158 MeV. We cannot rule out the possibility that a high density transition reappears at some temperature lower than the lowest temperature (129 MeV) that we have considered. Or perhaps the second critical endpoint goes away for light quarks. These possibilities we reserve for later investigation.
Finally, we take a look at an interesting observable introduced in [39], the ratio ξ/ξ 2 of correlation length and second momentum, which tests for the presence, in an effective Polyakov line action, of a spectrum of excitations contributing to two-point correlators. Our values are comparable to some of the results for pure SU(2) lattice gauge theory quoted in ref. [39], for more details and a thorough discussion see [31].

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 . This line has two critical endpoints, at T = 158 MeV and T = 220 MeV. It would also be interesting to see what happens to the second critical endpoint in a simulation with lighter quarks, or whether the expected transition line reappears, for m = 695 MeV, at temperatures below 129 MeV. We reserve these questions for later study.