A study of the sign problem for lattice QCD with chemical potential

We study the expectation value of the phase of the fermion determinant for Wilson lattice fermions with chemical potential. We use quenched SU(3) ensembles and implement a recently proposed exact dimensional reduction of the fermion determinant. Ensembles at several temperatures below and above the phase transition are studied and we analyze the role of the quark mass, the temperature, the volume and the topological sectors. We compare our numerical results to predictions from chiral perturbation theory.


Introductory remarks
With the increasing amount of experimental data on the QCD phase diagram, corresponding ab-initio lattice calculations become more and more important. However, when a chemical potential is introduced, lattice simulations face a serious challenge, the fermion sign problem. With non-zero chemical potential µ the fermion determinant det[D(µ)] is complex and cannot be directly used as a probability weight. Unless conceptually new ideas are developed, Monte Carlo simulations need to use various kinds of reweighting techniques.
The severeness of the sign problem, and thus the numerical effort for a reweighting strategy, may be characterized by the expectation value where e iθ is the phase of the fermion determinant det[D(µ)]. Having the extra factor of 2 in the exponent on the lhs. is convenient, since the expectation value of that phase may be written as a ratio of two determinants.
Recently the determinant phase (1) was addressed in several papers [1] - [6] using different analytical tools, such as chiral perturbation theory or random matrix models. The dependence of the sign expectation value e i2θ on the chemical potential, the volume, the temperature, the quark/pion mass and the topological sector was studied. On the lattice several individual results may be spotted [7] - [10], but a systematical analysis of the sign problem is still missing.
In this paper we attempt a small step towards a more complete analysis and study the determinant phase for quenched ensembles in a wide range of temperatures T and chemical potential values µ. Even in the quenched situation this is still a sizable task, but applying the recently proposed factorization formula [11,12] we are able to speed up the evaluation of the determinant phase considerably. We study the dependence of e i2θ on the parameters µ, T , the quark mass m, the volume and the topological sectors. In particular we also compare the behavior in the low-and high temperature phases of quenched QCD.

Technicalities
In our study we use quenched ensembles generated with the Lüscher-Weisz action [13] on lattices of size N 3 ×N T = 6 3 × 4, 8 3 × 4 and 10 3 × 4. The scale was determined in [14] based on the Sommer parameter. We generated ensembles for a wide range of temperatures between T = 210 MeV and T = 430 MeV, with the critical temperature for our action determined [15] to be T c = 300(3) MeV. Our statistics is between 500 and 2000 configurations for the smaller lattices and 100 configurations for the larger ones. An overview of our ensembles is given in Table 1.
For our gauge ensembles we determined the fermion determinant for the Wilson Dirac operator D(µ) with chemical potential µ (the lattice spacing is set to a = 1 here), We use the convention U −ν (x) = U ν (x −ν) † , introduce the abbreviationν for the shift vector in the ν direction, and κ is related to the bare mass m via κ = 1/(4 + m). The temporal boundary conditions for the fermions are anti-periodic, all other boundary conditions (gluons, spatial boundary conditions for the fermions) are periodic. In order to speed up the evaluation of the fermion determinant for many values of µ at the same time, we use the dimensional reduction formula developed in [11]. We here only very briefly sketch the idea of the construction and refer to [11,12] for the technical details: Applying a decomposition of the lattice into four domains, the fermion determinant may be rewritten in the form Here A 0 is a factor that is essentially a product of determinants for the terms of the Dirac operator restricted to the four domains. This factor is real and independent of the chemical potential µ. All of the µ-dependence comes from the second factor which has again the form of a determinant. However, the matrices H 0 , H ±1 live on only a single time slice and thus the evaluation of this second determinant is by a factor of N 3 T = 4 3 = 64 cheaper than the evaluation of the determinant in the original formulation. The matrices H 0 , H ±1 are made from products of propagators on the domains and are small enough such that they can be stored in memory. Then with (3) the determinant det[D(µ)] can be evaluated for several values of µ very efficiently. Actually, since A 0 cancels in (1), for the phase factor e i2θ only the dimensionally reduced determinant det[1 − H 0 − e µNT H +1 − e −µNT H −1 ] is needed. We typically use 16 to 26 different values for µ spaced with ∆µ = 0.05 (in lattice units).
For some of our ensembles we also evaluated the topological charge of the configurations. For that purpose the low lying eigenvalues of the overlap operator were computed and the topological charge was determined from the number of zero modes using the index theorem.

Qualitative dependence on mass and volume
We begin our presentation of the numerical results with a qualitative discussion of the behavior of e i2θ as a function of the quark mass and the volume. In the lhs. plot of Fig. 1 we show e i2θ for the  Similarly in the rhs. plot of Fig. 1 we compare the three volumes which we have at β G = 7.5 for m = 100 MeV. Here it is obvious, that the distribution is wider for small volume, i.e., increasing the volume makes the sign problem more severe.
We remark that although we here only show plots for ensembles in the low temperature phase, we confirmed that the same behavior is also found in the high temperature phase, i.e., small mass and large volume increase the sign problem. This subsection is only meant to qualitatively demonstrate the behavior of e i2θ as a function of the quark mass and the volume. However, at least in the low temperature phase one can go beyond that and compare the numerical results to quantitative analytical predictions. This will be done in the next section.
In this section we conclude with a check of the reliability of our numerical results. Since the phase of the determinant is a difficult to measure quantity -after all it is exponentially suppressed for large µ -such a consistency check is important. In principle the phase e i2θ is a complex number. However, due to the symmetry under time reflections e i2θ is real. Thus a check for the quality of the Monte Carlo result for e i2θ is to inspect the imaginary part of that quantity and to control whether it is compatible with zero. In Fig. 2 we show the imaginary part of the phase expectation value for two ensembles in the low-and the high temperature phases for three different quark masses. In all cases we find that the imaginary part is zero within error bars as expected.

Comparison to chiral perturbation theory
After the first round of a more qualitative assessment of e i2θ in the last subsection, we now focus on the low temperature regime where analytical results are available. Since for our ensembles we have m π L larger than 5 for all our ensembles 1 , we can compare to results obtained in the p-regime of chiral perturbation theory [1,2], while we cannot expect agreement with the random matrix theory results in the microscopic regime. The basis for our comparison are the chiral perturbation theory results for e i2θ discussed in [2]. In particular the result for the expectation value of the determinant phase in the quenched case reads e i2θ = e − g0(µ) + g0(0) , g 0 (µ) = V m 2 π π 2 β 2 ∞ n=1 1 n 2 K 2 (n m π β) cosh(n 2µβ) .
Here V is the 4-volume, β = 1/aN T is the inverse temperature, and K 2 denotes the modified Bessel function. The expression (4) is expected to describe the phase expectation value for µ < m π /2. In Fig. 3 we show again the results for the low temperature ensembles already used in Fig. 1, but now plot them as a function of µ in lattice units (which is more convenient for our fits). In the lhs. plot we performed individual one parameter (the pion mass m π ) fits for all three quark masses. The fit parameters for the three ensembles with bare quark masses m = 50, 100, 200 MeV, are am π = 0.812(25), 0.861(25) and 0.938(34), which corresponds to physical pion masses of m π = 682(21), 724(21) and 789(29) MeV 2 . Only data points with µ < m π /2, where (4) is expected to hold, were taken into account in the fit. It is obvious that the corresponding curves in the lhs. plot of Fig. 3 describe the data very well. Actually the values of χ 2 /d.o.f. are rather small (below 0.1) throughout.
The data in the rhs. plot of Fig. 3 allow for a much more stringent test of the chiral perturbation theory result (4). The three ensembles used there differ only by their volume, while the gauge coupling β G and the quark mass m were held fixed. Thus one expects that the three ensembles have the same pion mass (neglecting possible finite size effects). Since the volumes V are known for the three ensembles, one can attempt a single one parameter (the pion mass m π ) fit for all three volumes simultaneously. The result of this fit is shown by the full curves in the rhs. plot of Fig. 3 and obviously represents the data pretty well. The outcome for the fit parameter is am π = 0.729(52) which corresponds to a pion mass of m π = 675(49) MeV (note that here the gauge coupling is different from the one in the lhs. plot). We find χ 2 /d.o.f. = 0.485 which demonstrates that (4) describes the data for the different volumes very well. It is also interesting to note that the result for the pion mass from the fit (m π = 675 MeV) is rather close to what one would estimate from the quenched spectroscopy calculation [16] done for the Lüscher-Weisz action (m π = 660 MeV).

The role of the topological charge
An interesting effect has been discussed in [6]: In a random matrix calculation in the microscopic regime it was found that the phase expectation value e i2θ should depend on the topological charge Q of the gauge configurations. The distribution of e i2θ becomes wider as |Q| is increased, in other words, the sign problem is milder for higher charge sectors. The question which we address here is whether the topological effect is a specific feature of the microscopic regime or plays a role in general.
As already remarked, we determine the topological charge via the number of zero modes of the overlap Dirac operator using the index theorem. The overlap operator has a free parameter s which may be used to tune the locality [17]. This parameter, however, also influences the number of zero modes as it shifts the center of the overlap projection. We compare two different values s = 0.5 and s = 0.8 which give slightly different results for the topological charge. The agreement of the topological charges becomes better as β G is increased.
In Fig. 4 we show e i2θ as a function of the chemical potential in lattice units for our β G = 7.4, 8 3 × 4 ensemble, separating the gauge ensembles with respect to the topological charge Q. The lhs. plot is for s = 0.5, the rhs. for s = 0.8. The plots clearly indicate that there is no topological effect for e i2θ in our data. The curves are on top of each other within error bars, and the central values do not show the expected monotony (upwards trend with increasing |Q|). The same absence of the topological effect was also found for the other values of β G . One may speculate about the reason for the absence of the topological effect. It could well be that it is indeed seen only in the microscopic regime or only for a chirally symmetric Dirac operator such as the overlap, but of course also the analysis could be improved by using larger and finer lattices, where the topological charge is better defined. Since the cost for evaluating e i2θ rises tremendously with the volume, such a study must be left for the future.

Increasing the temperature
We now discuss the behavior of e i2θ as one increases the temperature. For the quenched theory the system undergoes a phase transition at T c ∼ 300 MeV, where the center symmetry is broken spontaneously. This breaking is, e.g., signaled by the Polyakov loop which vanishes below T c and has a finite value above T c . In the high temperature phase the system spontaneously selects one out of three possible center sectors, which are distinguished by the phase φ P of the Polyakov loop, φ P ∼ 0, 2π/3 or −2π/3. The fermion determinant on the other hand is not invariant under center transformations and becomes very small for the two complex sectors due to self averaging of the canonical determinants in the fugacity expansion [12]. Only for the real center sector the fermion determinant remains large, which in turn is the reason that in a dynamical simulation always the real center sector is selected. Consequently in our analysis we here consider the real center sector. We only take into account those gauge configurations where the phase φ P of the Polyakov loop obeys |φ P | < π/3. For the other sectors the fermion determinant becomes relatively small and the phase factor e i2θ is ill defined. In Fig. 5 we give an overview of the determinant phase for the 8 3 × 4 ensembles plotted as a function of the chemical potential in lattice units, comparing all our values of the gauge coupling. The critical value for the transition is roughly given by β G = 7.8. In Fig. 6 the same information is presented as a 3-d plot. The figures show clearly that qualitatively the distribution of e i2θ versus µ keeps the Gaussian type of shape as the temperature is increased. The distribution first seems to become more narrow with increasing β G , but from β G = 7.8 on widens again. However, it is important to keep in mind, that although in lattice units the volume remains fixed, the lattice spacing and thus the physical volume decrease with increasing β G . One expects that the shrinking physical volume has a mildening effect for the sign problem (see Subsections 3.1 and 3.2). The main message of Figs. 5 and 6 thus is of more qualitative nature: At the quenched transition there is only a gradual change in the distribution of e i2θ and we do not observe a dramatic qualitative effect for the sign problem.

Discussion
In this paper we have made a step towards a more systematical understanding of the fermion sign problem for lattice QCD with chemical potential. Using quenched ensembles the phase of the fermion determinant is analyzed for a wide range of temperatures below and above the phase transition and for several values of the chemical potential. We analyze the dependence of the determinant phase on the temperature, the chemical potential, the quark mass, as well as the topological charge and compare our results to chiral perturbation theory. We find that for all temperatures and values of the chemical potential the sign problem becomes harder with increasing volume and decreasing quark mass. In general the determinant phase has a Gaussian type of distribution as a function of the chemical potential for all temperatures we considered. At the deconfinement transition of the quenched theory we do not observe any dramatic qualitative effect for the sign problem. Concerning a possible dependence on the topological charge, we do not find such a topological effect in the regime we work at. However, here technical improvements, in particular larger and finer lattices, where the concept of topological charge is better defined, would be needed for a final answer on the fate of the topological effect outside the microscopic regime.