Polyakov line actions from SU(3) lattice gauge theory with dynamical fermions via relative weights

We extract an effective Polyakov line action from an underlying SU(3) lattice gauge theory with dynamical fermions via the relative weights method. The center-symmetry breaking terms in the effective theory are fit to a form suggested by effective action of heavy-dense quarks, and the effective action is solved at finite chemical potential by a mean field approach. We show results for a small sample of lattice couplings, lattice actions, and lattice extensions in the time direction. We find in some instances that the long-range couplings in the effective action are very important to the phase structure, and that these couplings are responsible for long-lived metastable states in the effective theory. Only one of these states corresponds to the underlying lattice gauge theory.


Introduction
Our approach to understanding the phase structure of QCD at finite densities is to map the theory onto a simpler theory, described by an effective Polyakov line action, and then to solve for the phase structure of that theory by whatever means may be available. At strong couplings and heavy quark masses the effective theory can be obtained by a strong-coupling/hopping parameter expansion, and such expansions have been carried out to rather high orders [1]. These methods do not seem appropriate for weaker couplings and light quark masses, and a numerical approach of some kind seems unavoidable. There are, of course, methods aimed directly at the lattice gauge theory, bypassing the effective theory. These include the Langevin equation [2] and Lefshetz thimbles [3]. In this article, however, we are concerned with deriving the effective Polyakov line action numerically, and solving the resulting theory at non-zero chemical potential by a mean field technique. In the past we have advocated a "relative weights" method [4,5], reviewed below, to obtain the effective theory, but thus far this method has only been applied to pure gauge theory, and to gauge theory with scalar matter fields.
Here we would like to report some first results for SU(3) lattice gauge theory coupled to dynamical staggered fermions. For an interesting alternative approach to determining the PLA by numerical means, so far applied to pure SU(3) gauge theory, see [6].

The Relative Weights Method
The effective Polyakov line action (henceforth "PLA") 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 x x,t 0) = 1), 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 Boltzman 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.
The relative weights method can furnish the following information about S P : Let U denote the space of all Polyakov line (i.e. SU(3) spin) configurations U x x x on the lattice volume. Consider any path through this configuration space U x x x (λ) parametrized by λ. Relative weights enables us to compute the derivative of the effective action S P along the path (dS P /dλ) λ=λ 0 at any point The strength of the method is that it can determine such derivatives along any path, at any point in configuration space, for any set of lattice couplings and quark masses where Monte Carlo simulations can be applied. The drawback is that it is not straightforward to go from derivatives of the action to the action itself, and in general one must assume some (in general non-local) form for the effective action, and use the relative weight results to determine the parameters that appear in that action.
In practice the procedure is as follows. Let denote two Polyakov line configurations that are nearby in U, with S L , S L the lattice actions with (1), where ... indicates that the expectation value is to be taken in the probability measure w.r.t. e S L . This expectation value is straightforward to compute numerically, and from the logarithm we determine ∆S P . Then (dS P /dλ) λ=λ 0 ≈ ∆S P /∆λ is the required derivative. The PLA inherits, from the underlying lattice gauge theory, a local symmetry under the transformation U x x x → g x x x U x x x g −1 x x x , which implies that the PLA can depend only on the eigenvalues of the Polyakov line holonomies U x x x . Let us define the term "Polyakov line" in an SU(N) theory to refer to the trace of the Polyakov line holonomy The SU(2) and SU(3) groups are special in the sense that P x x x contains enough information to determine the eigenvalues of U x x x providing, in the SU(3) case, that P x x x lies in a certain region of the complex plane.
Explicitly, if we denote the eigenvalues of U x x x as {e iθ 1 , e iθ 2 , e i(θ 1 +θ 2 ) }, then θ 1 , θ 2 are determined by separating (4) into its real and imaginary parts, and solving the resulting transcendental equations cos(θ 1 ) + cos(θ 2 ) + cos( In this sense the PLA for SU(2) and SU(3) lattice gauge theories is a function of only the Polyakov lines P x x x . We therefore compute the derivatives of the effective action, by the relative weights method, with respect to the Fourier ("momentum") components a k k k of the Polyakov line configurations The procedure is to run a standard Monte Carlo simulation, generate a configuration of Polyakov line holonomies U x x x , and compute the Polyakov lines P x x x . We then set the momentum mode a k k k = 0 in this configuration to zero, resulting in the modified configuration P x x x , where y y y P y y y e −ik k k·y y y e ik k k·x x x .
Then define where f is a constant close to one. We derive the eigenvalues of the corresponding holonomies U x and U x , whose traces are P x x x , P x x x respectively, by solving (5). The holonomies themselves can be taken to be diagonal matrices, without any loss of generality, thanks to the invariance under U x x x → g x x x U x x x g −1 x x x noted above. If we could take f = 1, then in creating P x x x , P x x x we are only modifying a single momentum mode of the Polyakov lines of a thermalized configuration. However, there are two problems with setting f = 1. The first is that at f = 1 and finite α there are usually some lattice sites where |P x x x |, |P x x x | > 1, which is not allowed. In SU(3) there is the further problem that at some sites the transcendental equations (5) have no solution for real angles θ 1 , θ 2 . So we are forced to choose f somewhat less than one; in practice we have used f = 0.8. The choice 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 , as described above, with respect to the real part a R k k k of the Fourier components a k k k .

A heavy-quark ansatz for S P
The problem is to derive S P from the derivatives ∂S P /∂a R k k k . Unfortunately there is no systematic procedure for doing this, and an ansatz for the effective action is required. For pure gauge theories we have assumed a bilinear effective action of the form This non-local coupling can be obtained from derivatives computed by the relative weights method. A test of the method and the ansatz (9) is to compare the Polyakov line correlator G(R) = P(x x x)P † (y y y) , where R = |x x x−y y y], computed in the effective theory with the corresponding correlator computed in the underlying lattice gauge theory. Excellent agreement was found in SU(2) and SU(3) pure gauge and gauge-Higgs theories [4,5,7]. Now we are interested in adding dynamical fermions, which break global center symmetry explicitly in the underlying lattice gauge theory, and the problem is to determine their contribution to the effective action. For heavy quarks the answer is known [8], and if we denote by S F the center symmetry-breaking piece of the effective action, then to leading order in the hopping parameter expansion, at non-zero chemical potential, we have where determinants can be expressed entirely in terms of Polyakov line operators, using the identities and where h = (2κ) N t , with κ the hopping parameter for Wilson fermions, or κ = 1/2m for staggered fermions, and N t is the extension of the lattice in the time direction. The power is p = 1 for four flavors of staggered fermions, and p = 2N f for N f flavors of Wilson fermions. It is possible to compute higher order terms in h in a combined strong-coupling/hopping parameter expansion [1], and of course fermion loops which do not wind around the periodic time direction will also contribute to the center symmetric part of the effective action.
Our proposal is to fit the relative weights data to an ansatz for S P based on the massive quark effective action, i.e.
where both the kernel K(x x x − y y 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 x x − y y y). We are aware that this is a strong assumption. There are two modest checks, however. First we can compare, at µ = 0, the Polyakov line correlators computed in the effective theory and the underlying gauge theory, and see how well they agree. Secondly, if it turns out that the h-parameter is very small even for quark masses which are fairly light in lattice units, then that is an indication that more complicated center symmetry-breaking terms are smaller still, and likely to be unimportant, at least at µ = 0. Finally, we do know qualitatively that an ansatz of this form satisfies the Pauli principle, in that the number density n of quarks per site will saturate, as µ → ∞, at the correct integer, which is n = 3 for three colors of staggered unrooted (p = 1) fermions. For these reasons we regard the ansatz (13) as a reasonable starting point for the relative weights approach, to be modified as necessary. Components of the wavevector k i = 2πm i /L are specified by a triplet of integer mode numbers (m 1 , m 2 , m 3 ), and in this work we have used triplets with lattice extension L = 16 in the three space directions. In calculating the center symmetrybreaking parameter h and momentum-space kernel K(k k k) at k k k = 0, we gain precision by carrying out the relative weights calculation at imaginary chemical potential µ/T = iθ. This is achieved by constructing U x x x ,U x x x as described above, and then making the replacements which are held fixed in the Monte Carlo simulation. The simulations are carried out for unrooted staggered fermions, corresponding to p = 1 in the heavy quark ansatz (13). The derivative of S P in (13) with respect to the real part a R 0 of the Polyakov line zero mode is then where If h 1, so that it is consistent to drop terms of O(h 2 ) and higher, then the derivative simplifies to The left hand side of this equation is computed numerically, for a variety of α, θ, by the relative weights technique. Plotting the data vs. α at θ = 0, we can find K(0) and h from the slope and intercept, respectively. However, a more accurate estimate of h is obtained by plotting the results vs. θ, at fixed α, and then extrapolating to α → 0. Having computed h and K(0), the next thing to do is to compute the kernel K(k k k) at k k k 0, and for this we can set the chemical potential to zero. We then have the derivative of the action with respect to non-zero modes a k k k of the Polyakov lines again dropping terms of order h 2 and higher. The left-hand side is computed via relative weights at a variety of α, and plotting those results vs. α, K(k k k) is determined from the slope.

Results for the PLA
In this initial study we have concentrated on parameters (β, quark mass ma, and inverse temperature N t in lattice units) which bring us close to, but not past, the deconfinement transition. In all cases we work on a 16 3 × N t lattice with staggered, unrooted fermions.

Wilson action, N t = 4
Figure 1a is a plot of K(k k k) vs. the lattice momentum for an underlying lattice gauge theory at β = 5.04, ma = 0.2. We have found in previous work [4], and find here also, that most of the data points can be fit by two straight lines The exception is one or two points at the lowest momentum, which do not fall on a straight line. If in fact K f it (k k k) would fit K(k k k) down to k L = 0, it would imply in position space that K(x x x −y y y) ∝ 1/|x x x −y y y| 4 . As in previous work, we interpret the deviation as implying a cutoff on the long range couplings, and define the position-space kernel with a long distance cutoff r max The cutoff r max is chosen so that, upon transforming this kernel back to momentum space, the resulting K(k) also fits the low-momentum data at low momenta. The result of this procedure is shown in Fig.  1b. The constant h = 0.033 is determined as explained in the previous section. The parameter h and the kernel K(x x x − y y y) are sufficient to specify the PLA, assuming the validity of the heavy-quark ansatz (13), and at zero chemical potential we may simulate both the PLA and the underlying lattice gauge theory (LGT) to compute and compare the Polyakov line correlators in each theory, shown in Fig. 2. We have also applied the relative weights method to the Lüscher-Weisz action, with the parameters listed above. Again most of the K(k k k) data points can be fit by two straight lines. However, there is a significant difference as compared to the previous cases at k k k = 0, where K(0) lies above, rather than below the straight line, as seen in Fig. 3a. Neglecting couplings between lattice sites beyond some separation r max will inevitably result in disagreement with the K(0) data point. In this case the strategy is to Fourier transform the two-line fit (20) to position space, with the modification that we identify K f it (0) with K(0), and dispense with a finite-distance cutoff at r max . The resulting kernel K(x x x − y y y) in the PLA couples each lattice site to every other lattice site. The result appears to be multiple metastable phases, which depend, in numerical simulations, on the initial configuration. In Fig. 3b we display our results for the Polyakov line correlator G(R) obtained from numerical simulations of the underlying lattice gauge theory and the Polyakov line action with a starting configuration initialized to P x x x = 0.3; and P x x x = 0.0; These results indicate that there are at least two phases in the PLA, confined and deconfined, which are stable over thousands of Monte Carlo sweeps. The Polyakov line correlator of the PLA in the confined phase is consistent with the correlator in the underlying lattice gauge theory, while the correlator in the deconfined phase is not. It seems that for the purpose of numerical simulations the effective action alone may be insufficient, and it may be necessary in some cases to supplement the PLA with a prescription for initialization of the SU(3) spin system.
The existence of multiple stable or metastable phases in the PLA is very clearly associated with the long-range couplings in the effective action. We have checked that if one simply truncates the range of bilinear couplings then the multiple phases disappear, and the result for the Polyakov line correlator is independent of the initialization. Of course, that arbitrary truncation also destroys the agreement of the correlators obtained in the PLA and the underlying lattice gauge theory.

Mean field solutions at finite density
We review here the mean field approach to solving the PLA at finite density, as explained in refs. [9] and [10]. The partition function corresponding to the action (13) is K(x x x − y y y)TrU x x x TrU y y y with D x x x (µ, TrU, TrU † ) and We then write Next introduce parameters u, v so that 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 u = TrU x , v = TrU † x . These conditions turn out to be equivalent to the stationarity of the mean field free energy. The leading mean field result is obtained by dropping E 0 , in which case the integrand of the partition function factorizes where A = J 0 v, B = J 0 u. The SU(3) group integral is known (see, e.g., [9]), where D s i j is the i, j-th component of a matrix of differential operators Putting everything together, with Z m f = exp[− f m f V /T ], the mean-field free energy/volume is In practice a computation of the mean field estimate f m f of the free energy requires a truncation of the sum over s in (32), an expansion in a 0 to finite order, and a check that the results are not sensitive to increasing the cutoff. We have found that restricting the sum over s to the range −3 ≤ s ≤ 3, and the expansion to a 0 to second order, is sufficient. The results for the examples we have considered in the last section, with the Wilson action and N t = 4, are qualitatively very much like the mean field results heavy quark cases, which were reported in [10]. The mean field solutions for TrU , TrU † and number density n for the case β = 5.04, ma = 0.2 are shown in Fig. 4.  The Lüscher-Weisz action at N t = 6, β = 7.0, ma = 0.3 is more interesting. There are multiple solutions of the mean-field equations (33), and the solution which is found by a search routine depends on the starting values for u and v. Initialization at u = v near zero gives the results shown in Fig. 5. Here there seem to be two clear phase transitions at finite density. If, however, the search routines begin at u = v = 1, then solutions correspond to the deconfined phase at µ = 0, and there is no transition found at any value of µ, as seen in Fig. 6. Ordinarily the stable phase corresponds to the phase with lowest free energy, and by this criterion (see Fig. 7) the solutions shown in Fig. 6 are selected. However, we have found that at µ = 0 this is not the phase which corresponds to the phase of the underlying lattice gauge theory. This of course raises the question of which metastable state corresponds to the state of the underlying gauge theory at finite density.

Conclusions
We have derived effective Polyakov line actions via the relative weights method for several cases of SU(3) lattice gauge theory with dynamical staggered fermions, and solved these theories at non-zero chemical potential by a mean field approach. At µ = 0 we find good agreement for the Polyakov line correlators computed in the effective theories and the underlying lattice gauge theories. We have also found, at µ = 0, that Polyakov line expectation values computed via mean field theory are in remarkably close agreement with the values obtained by numerical simulation, and this is probably due to the fact that each SU(3) spin is coupled to very many other spins in the effective theory, which favors the mean field approach. However, this non-local feature of the effective action also leads, in the most non-local case we have looked at (each spin coupled to all spins) to an unpleasant feature, namely, the existence of more than one metastable phase. These phases depend on the initialization chosen, and they persist throughout the numerical simulation, involving thousands of Monte Carlo sweeps. Since this is a phenomenon seen at µ = 0, it is not specifically tied to the sign problem, but rather to the non-locality of the effective action in certain cases. At µ 0 one must either find some criterion for choosing the phase which corresponds to lattice gauge theory, or else restrict the analysis of the Polyakov line action to cases where the couplings are comparatively short range. It should be emphasized that even if there are significant terms in the action which are ignored in the simple ansatz (13), and even if such terms were taken into account, there may still be multiple metastable phases if the bilinear couplings are long (or infinite) range. Whether, in such cases, some other simulation method (e.g. Langevin) could avoid the dependence of phase on initialization remains to be seen. It is also important to discover whether very long or infinite range couplings in the effective action are generic at small quark mass and small lattice spacings, or whether such cases are special and can be bypassed. See also [11][12][13].