Two-flavor lattice QCD with a finite density of heavy quarks: heavy-dense limit and “particle-hole” symmetry

We investigate the properties of the half-filling point in lattice QCD (LQCD), in particular the disappearance of the sign problem and the emergence of an apparent particle-hole symmetry, and try to understand where these properties come from by studying the heavy-dense fermion determinant and the corresponding strong-coupling partition function (which can be integrated analytically). We then add in a first step an effective Polyakov loop gauge action (which reproduces the leading terms in the character expansion of the Wilson gauge action) to the heavy-dense partition function and try to analyze how some of the properties of the half-filling point change when leaving the strong coupling limit. In a second step, we take also the leading nearest-neighbor fermion hopping terms into account (including gauge interactions in the fundamental representation) and mention how the method could be improved further to incorporate the full set of nearest-neighbor fermion hoppings. Using our mean-field method, we also obtain an approximate (μ, T) phase diagram for heavy-dense LQCD at finite inverse gauge coupling β. Finally, we propose a simple criterion to identify the chemical potential beyond which lattice artifacts become dominant.


Introduction
Lattice QCD at non-zero baryon density suffers from the so-called "sign problem" [1]. It also displays a new type of lattice artifact: because of the Pauli exclusion principle, the lattice cannot support more than one quark per site (and per spin, color and flavor). The main goal of this paper is to characterize the effects of this new discretization error, so as to better disentangle genuine physics from lattice artifacts.
We first briefly recall the foundations of lattice QCD in the presence of a chemical potential µ and some well-known features of the sign problem which appears when µ = 0.

JHEP02(2016)051
We then focus on the half-filling point (section 2) and the associated fermion saturation. Using the heavy-dense limit (section 3), we elucidate the properties of the half-filling regime and analyze, using a mean-field calculation which takes leading nearest-neighbor fermion hoppings and gauge interactions into account (appendix B), how far some of these findings generalize to full QCD. Finally, using our mean-field method, we obtain an approximate (µ, T ) phase diagram for heavy-dense LQCD, which can be compared with current state of the art complex Langevin investigations [3].

QCD at finite density
In this paper, we adopt the Wilson discretization of the Dirac operator. For a single flavor at chemical potential µ it reads: where (x, y) are multi-indices labeling different lattice sites, (I, J), (a, b) are color and Dirac indices respectively and U ν,IJ (x) is the (I, J) component of a link variable that points from site x to (x + ν). The hopping parameter κ is related to the bare lattice quark mass m by κ = (2 m + 8) −1 , and the boundary conditions are assumed to be periodic in spatial directions and anti-periodic in the time direction.
Using matrix notation, we write (1.1) simply as D(µ), keeping its dependency on the gauge field implicit. From the definition (1.1), it follows that where µ * is the complex conjugate of µ, which shows that which can be used, together with Wilson's lattice gauge action, (1.7) to sample gauge configurations by Monte Carlo. This generalizes to the case of N f flavors as long as the degeneracy factor for each quark mass is even, i.e. as long as there is always an even number of flavors that have the same quark mass.

The sign problem
For Re(µ) = 0, the fermion determinant (1.3) becomes in general complex and loses its probabilistic interpretation. In order to continue sampling the gauge field by Monte Carlo, one is then forced to invoke reweighting techniques. For example, if the theory contains an even number N f of mass degenerate flavors, one can use |Det(D(µ))| N f = Det(D(µ)) N f /2 Det(D(−µ)) N f /2 instead of Det(D(µ)) N f to sample the gauge field and treat the complex phase, . (1.10) Although the reweighting procedure (1.9) is in principle exact, its applicability is limited by the statistical fluctuations of the phase R[U ]. To see this, consider the expectation value of the phase, where ∆f = f − f q is the difference in the free energy densities of the unquenched and the quenched theory and L and N t are respectively the spatial and temporal system sizes.

JHEP02(2016)051
For sufficiently large systems, ∆f becomes approximately independent of the system size and (1.11) shows that R q then decays exponentially with increasing system size. As per definition the modulus of R[U ] remains equal to 1, this shows that the fluctuations in the reweighted observable, O[U ] R[U ] / R q , grow exponentially with the system size and so does the statistical error. The dependency of the latter on the number of measurements is known to be proportional to 1 / √ # meas.. This means that the statistics required to obtain equally accurate results for different system sizes must scale like required statistics ∝ e 2 L 3 Nt ∆f . (1.12) The exponential growth of the required statistics as a function of the system size in (1.12) is a manifestation of the so-called sign problem which is generic in the numerical treatment of strongly correlated fermion systems at finite density.

Half-filling and saturation
As mentioned above, the sign problem is generic in the numerical study of strongly correlated fermion systems at finite density and appears therefore also in solid state physics when working with e.g. the fermionic Hubbard model. There, finite density simulations are often carried out at the so-called half-filling point, i.e. at the value of the chemical potential where the fermion density assumes half its maximum value. Due to the particle -hole symmetry that exists at this point, the sign problem disappears and it is possible to simulate much larger systems than at generic non-zero values of the chemical potential. The half-filling point also exists in lattice QCD (LQCD). The reason for this is of course that the quarks, as they are fermions, are subject to the Pauli exclusion principle, which implies that on each site, each of the 2 × N f × N c fermionic states can be occupied only once. For two flavor QCD (N f = 2, N c = 3), the maximal number of fermions per site is 12 and the half-filling point should therefore correspond to a fermion number density of 6. Figure 1 illustrates this for the case of a tho flavor system with an isospin chemical potential (i.e. a phase quenched two flavor system): it shows the isospin number density as a function of the chemical potential µ for different system sizes. The figure also illustrates that the location of the half-filling point is independent of the system size as the different curves (corresponding to different system sizes) all cross in precisely the same point at n I = 6, µ ≈ 1.3. In figure 2 we also show the condensates ψ ψ and ψγ 4 τ 3 ψ as a function of isospin-µ for the same set of system sizes. Also here we see volume independency at half-filling as the different curves cross in a single point at µ ≈ 1.3, but in contrast to the situation with the isospin density, the values of the condensates at half-filling are not half-way between their corresponding values at µ = 0 and at saturation µ → ∞. To see where the volume independency comes from, we define the local isospin density,  Figure 2. Quark (left) and "isospin" (right) condensate as a function of isospin µ for different system sizes, recorded at κ = 0.15 and β = 5.0. All curves cross in a single point, the half-filling point at µ ≈ 1.3, but in contrast to the situation with the isospin density in figure 1, the values of the condensates at half-filling are not in the middle between their values at µ = 0 and µ → ∞.
where tr c,d is the trace with respect to color and Dirac indices, such that and consider its average space-time variance, given by This quantity is shown on the left-hand side of figure 3 for the same ensembles that were used in figure 1. As can be seen, (2.4) is essentially zero at the half-filling point, which means that the half-filling is realized very homogeneously. As, according to the right-hand  Figure 3. Average isospin density space-time variance eq. (2.4) (left) and isospin density ensemble variance (right) as functions of µ for different system sizes, recorded at κ = 0.15 and β = 5.0. The figures show that at the half-filling point at µ ≈ 1.3, both the space-time variance and also the ensemble variance of the isospin density are essentially zero, which suggests that the isospin density becomes independent of the gauge fields for this particular value of µ. side of figure 3, also the ensemble variance of the total isospin density (2.3) is approximately zero at the half-filling point, we can conclude that the homogeneity of the local isospin density (2.2) at this point is due to the fact that its dependency on the gauge field is highly suppressed at half-filling, which explains the volume independency.
In figure 4 we also show the average sign (or average phase) (1.11) for a two-flavor system as a function of the chemical potential (left) and as a function of the isospin number density (right), for two different sets of simulation parameters: once for a system with L 3 N t = 4 4 , inverse coupling β = 5.0 and hopping parameter κ = 0.15 and once for one with L 3 N t = 5 4 , β = 5.3 and κ = 0.175. As can be seen, for both systems the average sign has a local maximum at the half-filling point, where it almost reaches a value of one again. However, for the larger system a clear deviation from one is visible. The right-hand part of figure 4 also shows an approximate symmetry about the half-filling point.
Although the average sign is not exactly one at half-filling, the fact that it has a local maximum, implies that isospin number density and quark number density must be the same at this point, as we have , (2.6) which vanishes at an extremum of R q (µ). This is of course also true at the minima of R q (µ), which means, that although the sign problem is worst at these points, the quark number density could there be obtained with high accuracy. The problem is, however, to determine the location of these points, as they will not be volume independent as is the case for the half-filling point. In figure 5 Figure 4. Average sign as a function of µ (left) and as a function of the isospin number density (right) for two different sets of hopping parameter κ, inverse coupling β and system size. The figures show that the average sign has a maximum at half-filling where its value is close to one. For the larger system there is however a clearly visible deviation from unity.  . Difference between average quark number density and average isospin number density as a function of µ for the same ensembles used in figure 4. The quark number density was obtained by reweighting.  The dotted lines correspond to fits of the form n I,sat − n I (µ) = a e Nt Q (µ) to the data for µ > 1.5μ, whereμ = − log(2κ) is the half-filling value of µ in each ensemble. As can be seen, both fits yield Q ≈ −2, although the simulation parameters are quite different. figure 4. The density difference was obtained by reweighting the quark number density from the quenched to the un-quenched system and then subtracting the isospin number density.
We mentioned that, in the Hubbard model, the sign problem disappears at half-filling due to a particle-hole symmetry. We could now ask if also in LQCD such a particle-hole symmetry is the origin of the (almost) disappearing sign problem at half-filling. A first indication that the answer to this question could be yes, comes from figure 6, which shows, again for a two-flavor system, coupled to an isospin chemical potential µ, the quantity n I,sat − n I (µ) for two different ensembles, where n I,sat = lim µ→∞ n I (µ) is the saturation value of the average isospin density n I (µ) at large µ. We could think of n I,sat − n I (µ) as being the density of "holes in the saturated state".

JHEP02(2016)051
For the isospin density itself we would expect n I (µ) ∝ e Nt Q µ for µ μ, wherẽ µ = − log(2κ) is the half-filling value of the chemical potential and Q = 2 is the charge of an elementary excitation (charged meson). The fits in figure 6 now show that n I,sat − n I (µ) behaves in precisely the same way for µ μ, but with Q = −2, which means that the "hole-excitations" from the saturating state at large chemical potential behave just like the anti-particles of the elementary excitations of the theory at small chemical potential. The strong interaction treats therefore fermion holes in the saturated state in a similar way as the fermions themselves in the vacuum, such that the holes also have to form gauge invariant combinations of the form of mesons and baryons.
Having this symmetry between particles at zero density and holes at saturation, makes it plausible that at half-filling, where particles and holes appear with equal densities, they indeed behave symmetrically. The same should be true for the non-phase-quenched system where the elementary excitations would be baryons with Q = 3 at zero density and holes of the form of anti-baryons with Q = −3 at saturation. Again this suggests that at half-filling, where baryons and baryon-holes have equal density, we should have particle-hole symmetry.

Dominance of lattice artifacts: a simple criterion
In contrast to the systems studied in solid state physics, where the lattice is physical and half-filling describes a physical state as in the fermionic Hubbard model, the lattice in LQCD is jut a regulator which has to be removed (e.g. in numerical works, by doing continuum extrapolation) in order to extract real physics. The half-filling state of a LQCD system is therefore just an unphysical lattice artifact and should be of no relevance for the continuum physics.
As long as we are running simulations sufficiently far away from the continuum limit (such that the maximal fermion density on the lattice (in physical units) is far below the density where the Pauli exclusion principle would become important in the continuum), the first appearance of a minimum in the average sign as a function of increasing chemical potential (or as a function of increasing isospin number density) as shown in figure 4, indicates the point where lattice artifacts start to become dominant. This becomes clear by noting that the average sign is a measure for the overlap between a system with a quark chemical potential and the corresponding phase-quenched system (with an isospin instead of a quark chemical potential). As the right-hand side of figure 4 shows, the average sign drops dramatically as soon as the isospin density becomes non-zero, reflecting the fact that a pion condensate develops in the isospin system, leading to a ground state which is rather different from that of the corresponding system with a quark chemical potential. The fact that the overlap between the two systems starts to increase again with increasing µ (i.e. that the average sign starts to increase again) indicates, that the influence of the Pauli exclusion principle starts to become dominant on the lattice: the fact that the quark chemical potential favors the u and d quarks, while the isospin chemical potential favors u and d, starts to become less important as for example the "d-hole" in the system with the quark chemical potential starts to play the role of the d in the isospin system and vice versa. This happens at densities significantly below half-filling, and measurements of observables taken at larger densities should be considered with caution.

JHEP02(2016)051 3 Heavy-dense analysis
In the following section we will analyze the half-filling point of LQCD in the heavy-dense limit and at strong coupling. These simplifications will allow us to understand the origin of the properties of the half-filling point.

The heavy-dense approximation
The goal of this section is to obtain an expression for the Partition function in the so-called heavy-dense approximation at strong coupling, which can be integrated analytically. The heavy-dense approximation usually consists of taking the limit κ → 0, µ → ∞, keeping κ e µ finite while terms proportional to just κ or κ e −µ vanish. In our case, it would be more appropriate to call it the static quark limit, as we will just drop the spatial hopping terms from the Dirac operator (1.1), but not the one corresponding to backward hopping in time, which is proportional to κ e −µ . The resulting fermion determinant is well-known as the leading term in the spatial hopping expansion that was already used in [6] to derive an effective model for QCD or in [2] to obtain a holomorphic action for complex Langevin simulations of LQCD.
Let us start by expanding the determinant of the Dirac operator (1.1) in terms of closed fermion loops C l 0 of length l 0 [4]: with M C l 0 = H x 1 x 2 · · · H x l 0 x 1 being the matrix product of hopping terms H xy = (S + T ) xy along the contour C l 0 of length l 0 , where S and T are the spatial and temporal hopping matrices defined in (1.1). The subscripts "c" and "d" of "det" and "tr" indicate that the operators act on color and Dirac space only. To get from the second to the third line in (3.1), we have used that k can be written as k = l 0 n, where l 0 is the length of the contour C l 0 , and n the number of windings around C l 0 . The change in the denominator (k → n instead of k → l 0 n) comes from the fact that Tr H l 0 n contains tr c,d M C l 0 n exactly l 0 times, once for every possible base point (or starting point) of a loop on C l 0 , which cancels the l 0 in the denominator. As already mentioned, the static quark approximation now consists of dropping the spatial hopping matrix S from (1.1) and (3.1), which ultimately means that M C l 0 is nonzero only if C l 0 wraps around the temporal direction, i.e. l 0 = n t (we use from now on JHEP02(2016)051 a lower-case n t to represent the temporal system size) and M Cn t has to be one of the n x · n y · n z matrices − e µ nt P x ⊗ (1 − γ 4 ) nt = T x,x+ 4 · · · T x+(nt−1) 4,x (3.2) or their time-reversed versions where P x is the Polyakov loop at spatial position x. The minus signs on the left-hand sides of (3.2) and (3.3) come from the anti-periodic boundary conditions in temporal direction. With this approximation, the determinant (3.1) reduces to a product of single site determinants: More precisely, the last line of (3.4) shows that for each site, we get four factors, each in the form of a determinant in color space: the first two correspond to the two distinct states of a spin-1/2 particle while the second two are due to its anti-particle. These color-determinants can be re-written in terms of traces of powers of the Polyakov loop: By using that the characteristic equation for a matrix X ∈ SU(3) reads and, according to the Cayley-Hamilton theorem, X satisfies its own characteristic equation, χ X (X) = 0 (in the spectral sense), it follows that (3.7) and therefore tr c X 3 = 3 + tr c (X) tr c X 2 − tr c X † tr c (X) , (3.9) tr c X † = tr 2 c (X) − tr c X 2 /2, (3.10) which can be used to simplify (3.5) further: Finally, expressing the traces in (3.11) in terms of the eigenvalues of the Polyakov loop P x , we arrive at det 2 and similarly, 14) The single site fermion determinants in (3.4) reduce therefore to the simple form: Using the Haar measure for the trace of the Polyakov loop, (3.17) and neglecting the plaquette action (i.e. considering the strong coupling limit β = 0), we can write down a heavy-dense partition function for QCD with two flavors, u, d:

JHEP02(2016)051
where κ u/d and µ u/d are the hopping parameter and the chemical potential corresponding to the u/d flavor, Λ is the spatial lattice with |Λ| sites, and Z s (µ u , µ d , κ u , κ d , n t ) is the single site partition function.
The integral (3.18) can be evaluated explicitly but results in a rather lengthy expression. For the purpose of illustration its explicit form is shown in appendix A for the case of degenerate quark masses, κ u = κ d = κ, once with an isospin (µ = µ u = −µ d ) and once with a quark (µ = µ u = µ d ) chemical potential.

Heavy-dense staggered fermions
The Staggered fermion operator is given by with the Staggered phase η µ (x) = (−1) 3 x 1 +2 x 2 +x 3 and the bare lattice fermion mass m.
The full static quark fermion determinant for (3.19) is more difficult to obtain than the Wilson fermion analogue (3.4). The reason is that Staggered fermions allow for retracting fermion loops, such that even without spatial hoppings, several distinct temporal loops are possible. However, if we are just interested in the heavy-dense limit, where µ and m are assumed to be large while e µ−log(2 m) is of order one, the leading term reads (3.20) Comparing (3.20) with (3.4) reveals that the Staggered heavy-dense determinant (3.20) can be obtained from the heavy-dense limit of (3.4) by replacing 2 κ by 1/(2 m) and dropping the square of the color-determinant.

Analytic expressions for some observables
In this section we derive some observables with respect to the system described by (3.18).

Average sign
In terms of (3.18), with κ u = κ d = κ, the average sign is computed as where φ is the complex phase of the single flavor determinant (3.4), and n x n y n z is the spatial system size.

Isospin and quark number density
The isospin density is obtained from (3.18) by

Heavy-dense quark propagator
The heavy-dense quark propagator is most easily obtained by using that [6] and We then obtain: where P y (y 4 ) is the Polyakov loop with base point y and P y (x 4 , y 4 ) is a Polyakov line at spatial position y that, for x 4 = y 4 , starts and ends at euclidean times x 4 , y 4 respectively (by going in positive time direction), while P y (x 4 , x 4 ) = 1. The inverse of the matrix 1 + (2 κ e µ ) nt P y (y 4 ) can easily be carried out (e.g. with Cayley-Hamilton and Leverrier-Faddeev), provided the matrix is not singular: where the determinant in the denominator is given by (3.11). Similarly we get: and therefore with (3.25): The heavy-dense quark propagator (3.32) could be used to get analytic expressions for average meson and baryon propagators from which one could determine the corresponding masses. But due to the projection properties of (3.32) in Dirac space, the mass spectrum will be highly degenerate. We will therefore use (3.32) in the following just to facilitate the derivation of some other observables.

Local isospin number density
The local isospin number density in the heavy-dense case can be obtained, by starting from equation (2.2) and replacing the full Dirac operators and propagators by their corresponding heavy-dense counterparts, i.e. just dropping the spatial hopping terms from (1.1) and using (3.32) as the corresponding propagator. We then obtain: where after the second equality sign it was used that (1 ± γ 4 ) are orthogonal projectors and that for example i.e. the link variable U 4 (x) coming from the term in the Dirac operator that is not projected out by the (1 − γ 4 ) of the propagator, closes the Polyakov line P x (x 4 + 1, x 4 ) to a Polyakov loop P x (x 4 ), starting and ending at x. Analogously:

Condensates
The quark condensate for a two-flavor system, is in the static-quark limit obtained by taking the expectation value of x is the static-quark propagator (3.32) for y = x and κ u , κ d and µ u , µ d are the hopping parameters and chemical potentials for the two flavors.
Analogously the isospin condensate, is obtained as the expectation value of where in the first line 1 c,t is the identity with respect to color and time indices.

Half-filling point
Having obtained the heavy-dense partition function for two-flavor LQCD in the strong coupling limit, and defined some basic observables, we now use these results to study the origin of the properties of the half-filling point in the heavy-dense/strong coupling limit, and check, how far our findings could also apply to full LQCD.

Location of the half-filling point
To find the location of the half-filling point, we consider the isospin density (3.22) in its integral form: For sufficiently large values of µ, the second term within the curly brackets in (3.40) can be neglected; setting µ = − log(2κ), this term is of order O (2κ) 2nt while the first term is of order 1 and becomes completely independent of P , as and we obtain for the isospin density (3.40): which means that µ =μ ≡ − log(2κ) corresponds to the half-filling point, independently of n t (up to O (2κ) nt corrections), as is visualized in figure 7, where the isospin number density (3.22) is shown as a function of µ for different system sizes and it can be seen that the curves, corresponding to different system sizes, all cross at the half-filling point, just as in the case of full LQCD shown above in figure 1. From (3.40) and (3.41) it also follows that the dependency of the local isospin density (3.33) on the gauge field is highly suppressed at half-filling: gauge field dependent terms are suppressed by a factor of O (2κ) 2 nt such that the space time variance of (3.33) nearly vanishes at half-filling, again just as in the case of full LQCD shown above in figure 3. Nevertheless, in full LQCD, the half-filling point is shifted towards a slightly larger value of µ: while in the heavy-dense case, we haveμ = − log(2 κ) ≈ 1.2 for κ = 0.15, we find in full LQCD, for the same value of κ and with β = 5.0, that the half-filling point is located at µ ≈ 1.3. This shift is due to spatial fermion hopping which causes corrections to the quark mass. Corrections coming from spatial hoppings and a finite value of β will of course also appear in the expressions for observables and give rise to deviations that could be much larger than O (2 κ) 2 nt . These corrections are considered in sections 3.5 and 3.6.

Particle-hole symmetry and average sign
Consider the single flavor fermion determinant (3.15) under the transformation µ → 2μ−µ that corresponds to reflection about the half-filling point: where the last line of (3.43) is true for µ ∈ [0, 2μ]. Comparing (3.43) with Det(D(µ, κ, n t )) * = Det(D(−µ, κ, n t )) = det 2 we see that they agree up to an irrelevant gauge field independent pre-factor (which is unity at the half-filling point) and terms which are suppressed by at least a factor of (2κ) nt . As the complex conjugate of the fermion determinant can be interpreted as resulting from a charge conjugation of the original fermion fields, this shows that the half-filling state possesses an exact T = 0 particle-hole symmetry also in the heavy-dense/strong coupling limit of LQCD. We can also check the behavior of the heavy-dense quark propagator (3.32), under the transformation µ → 2μ − µ. For y 4 = x 4 and µ ∈ [0, 2μ] we find: which by noting that and remembering that shows that: i.e. up to a minus sign and terms of higher order in (2κ), we find that for y 4 = x 4 , reflecting µ aboutμ turns the quark propagator into its Hermitian conjugate. For the equal time part, y 4 = x 4 we find in a similar way: For the two condensates introduced in section 3.2.5, we then find for the phase quenched case: and at half-filling exactly in the middle between their corresponding values at µ = 0 and for µ → ∞.
In the un-quenched case, the relations would be: and Returning to the single flavor fermion determinant (3.15), we can see from (3.43) and (3.44) that it is essentially real at µ =μ and the average sign (3.21) should therefore be 1 + O (2κ) nt . This can also be verified more directly by using (3.15) to show that for µ =μ, we have: det c 1 + (2κ eμ) nt P = 1 + (2κ eμ) nt tr c P + (2κ eμ) 2 nt tr c P † + (2κ eμ) 3 nt = 2 + tr c (P ) + tr c (P † ) = 2 (1 + Re(tr c (P ))) , (3.58) i.e. we find again that the dominant part of the single flavor fermion determinant (3.5) at half-filling, det 2 is real and positive up to terms of order O (2κ) 2 nt . The average sign (3.21) in the heavy-dense limit is shown in figure 10 as a function of the chemical potential (left) and as a function of the isospin number density (right) for the same system sizes and values of the hopping parameter κ that were used to generate the corresponding plots in full LQCD shown above in figure 4. By comparing the two figures, it can be seen, that in the heavy-dense case, the average sign is much more symmetric about the half-filling point and at the half-filling point itself, it deviates much less from being unity than in full LQCD, which is again due to the absence of spatial fermion hoppings, such that the deviation is really just of order O (2κ) 2 nt . 3.4 Mean-field method for heavy-dense QCD at finite gauge coupling β So far we have only considered heavy dense LQCD in the strong coupling limit, β = 0.
In this section, we would now like to get a glimpse of finite β effects by introducing an effective nearest-neighbor Polyakov loop action, which is motivated by the leading terms of the character expanded gauge field Boltzmann factor, as derived for example in [6], and then using a mean-field approximation. Our mean-field treatment is similar to one of the approaches considered in [8] but differs from the mean-field calculation in [7] as we are considering the Polyakov loops as our effective degrees of freedom while in [7], this role was taken by spatial links. The nearest-neighbor Polyakov loop action is obtained as follows: we start with a character expansion of the SU(3) Yang-Mills Boltzmann factor (cf. [5], chapter 3.4), where the product is over all plaquettes p, χ r (U p ) is the character of the SU(3) group element U p (corresponding to the plaquette p) in the representation r, and a r (β) = c r (β)/c 0 (β) with c 0 (β) being the expansion coefficient of the trivial character χ 0 (U p ) = 1.
We now drop the summation over r and keep just the terms corresponding to the fundamental representation. After integrating out the spatial links, we are then left with [6]: where the products run over all pairs of nearest-neighboring sites, i, j . If we now require that the Boltzmann factor corresponding to the desired effective nearest-neighbor Polyakov loop action should, to lowest order, be proportional to (3.61), we find that the effective JHEP02(2016)051 action should be given by: where L i = tr c (P i ). At this point it should be mentioned that in contrast to the situation with the U(1) or SU(2) Yang-Mills action, where the coefficients of the character expansion can be written in closed form in terms of Bessel functions, no closed form is known for these coefficients in the case of SU (3). To compute a f (β), we made use of the power series representations for c 0 (β) and c f (β), given in [5]  We now proceed by applying the mean-field approximation to (3.62), i.e. we write L i = L + δL i , keep only terms of order O δL , and then write again δL i = L i − L (and proceed analogously for L * i ), which, after dropping constant terms, leads to the following single site mean-field action: which we use to define an additional probability weight, w L, L * , L, β, n t = e 6 a n t f (β) (L * L + L * L) , (3.66) that has to be included in the single site partition function defined in (3.18), i.e.: For general values of µ u , µ d , κ u and κ d , the product of the two fermion determinants in (3.67) is usually complex. As a consequence L and L * differ and it is a subtle issue how to proceed in a mean-field treatment [8]. To bypass the subtleties associated with L * = L when setting L = L and L * = L * , we define L to be the mean-field value of L with respect to the phase quenched system, where L * q = L q holds and the identification L = L q , L * = L * q therefore leads to L * = L. The mean-field L is therefore determined by solving the following self-consistency equation: |Det(D(θ 1 , θ 2 ; µ u , κ u , n t )) Det(D(θ 1 , θ 2 ; µ d , κ d , n t ))| w L, L(θ 1 , θ 2 ), β, n t , (3.68) where w L, L(θ 1 , θ 2 ), β, n t ≡ w L, L, L(θ 1 , θ 2 ), β, n t and where is the phase-quenched partition function. After one has determined the stationary L with respect to (3.68) for a particular set (µ u , µ d , κ u , κ d , n t ) of parameters, one can reweight to the non-phase-quenched system by using L µ u , µ d , κ u , κ d , n t in (3.67) when computing expectation values of observables.
This difference between L * µ u , µ d , κ u , κ d , n t and L µ u , µ d , κ u , κ d , n t can clearly be seen in figure 11 where we show the two quantities for the mass degenerate case (i.e. κ u = κ d ) as a function of a quark chemical potential µ = µ u = µ d .
It has been mentioned in [8] that in a mean-field treatment, the reweighting method can be no more than an approximation scheme: while in a full Monte Carlo simulation, reweighting is exact in the limit of infinitely many sampled configurations, a mean-field treatment relies only on the most dominant configuration. In our opinion this argument does not hold as in a mean-field treatment, the expectation value of an observable is determined on the active site, usually by integrating exactly over all possible values of its configuration variables, and reweighting should therefore work perfectly well! Setting L * = L, the integrals in (3.67) and (3.69) can be carried out exactly, leading to solutions in terms of infinite sums of modified Bessel functions of the first kind: respectively, where the coefficients C l 1 ,l 2 µ u , µ d , κ u , κ d , n t are given by Unfortunately, the evaluation of (3.72) and (3.73) to the required accuracy is numerically rather expensive and it turned out to be more efficient to use direct numerical integration to solve for L in (3.68) and for computing observables. Expressions (3.72) and (3.73) were therefore merely used for cross-check purposes. An expression for (3.72) in the more general case of L * = L could be found along the lines of [9] where the mean-field action (including correction terms) is derived for the case of a pure SU(N ) Polyakov line model with a chemical potential.

Mean-field heavy-dense phase diagram and finite β effects
In figure 12 we show the average Polyakov loop for the mass degenerate two-flavor case as a function of isospin chemical potential µ = µ u = −µ d and inverse coupling β for five different temperatures, n −1 t , as obtained with our mean-field method described in the previous section. As can be seen, for sufficiently large β a deconfinement transition occurs and the value of β at which this transition happens becomes smaller with increasing temperature (decreasing n t ).
A finite value of the inverse gauge coupling β has also an effect on the average sign, which is illustrated in figure 13. There we show mean-field results for the average sign as a function of µ = µ u = µ d for a mass degenerate (κ = 0.15) heavy-dense two-flavor system of spatial size V = L 3 = 4 3 and temporal extent N t = 2, 3, 4, 8. The left-hand part shows the strong-coupling limit, β = 0, while the right-hand part shows the β = 5.0 case. For static quarks, the location of the half-filling point is not affected by the finite β value, but the overall temperature dependency of the average sign clearly changes: at finite β the sign-problem becomes weaker with increasing temperature and almost disappears after the deconfinement transition.
Note that the computation of the average sign with our mean-field method comes with a slight complication: the active site couples to its six nearest neighbors and there are therefore 6 interaction terms entering the mean-field action. But the ratio of links to sites should be 3 : 1 on a periodic 3-dimensional lattice, which means that by using simply a formula of the form of (3.21) to compute the average sign with our mean-field partition functions (3.67) and (3.69), i.e. e 2i φ (µ, κ, β, n t , n x n y n z ) = Z s µ, κ, n t , L(µ, κ, n t ) Z s,q µ, κ, n t , L(µ, κ, n t ) smallest representative subset of a larger system, in which for example even sites are active and odd sites are passive or vice versa, and the computed sign will correspond to such a system). This factor is obtained by taking the ratio of the average fermion determinants in the non-phase-quenched and phase-quenched case. The average sign then becomes: e 2i φ (µ, κ, β, n t , n x n y n z ) = Z s µ, κ, n t , L(µ, κ, n t ) Z s,q µ, κ, n t , L(µ, κ, n t ) Det(D) Det(D) µ, κ, n t , L(µ, κ, n t ) |Det(D) Det(D)| q µ, κ, n t , L(µ, κ, n t ) nx ny nz /2 , (3.76) where Det(D) is the static quark determinant (3.4) and the factor of 1/2 in the exponent in (3.76) is there because the term inside the bracket is now the reweighting factor for two sites. Figure 14 finally shows the average Polyakov loop as a function of temperature and isospin (left) or quark (right) chemical potential i.e. it essentially shows the (T, µ) phase diagram of our simplified model in the phase-quenched and un-quenched case: one can read off how the pseudo-critical temperature for the deconfinement transition changes as a function of the chemical potential. The right-hand part of figure 14, corresponding to the un-quenched case, can be compared, e.g. to the phase diagram in [3, figure 2], obtained by complex Langevin. In order to simplify the comparison, figure 15 shows the data from figure 14 using the same scales used in [3, figure 2], assuming that the lattice spacing of a = 0.15 fm, determined in [3], should apply equally well to our system, as we used the same simulation parameters κ = 0.04, β = 5.8. In our effective model, we take the gauge field only in the fundamental representation into account. But at high temperature, effects coming from higher representations are much less suppressed than at low temperature (the coefficients a r (β), if included in the effective Polyakov loop action, would appear only with small exponents, leading to a weaker suppression of the higher order terms), which explains why in our figure the deconfinement transition is generally shifted towards larger temperatures compared to [3, figure 2].

From the heavy-dense to the full fermion determinant
The heavy-dense fermion determinant is the leading term in the so-called spatial hopping expansion of the full fermion determinant which, for one flavor, can be obtained as follows: where S and T are the spatial and temporal hopping terms as defined in (1.1) and we have added a subscript s to the spatial hopping parameter, just for book-keeping purposes. The matrix (1 − κ T ) −1 is the heavy dense quark propagator given by (3.32). We can now proceed in a similar way as in (3.1) to write (3.77) in terms of a product of smaller JHEP02(2016)051 determinants of "loop matrices", but this time, the loops are purely spatial and the matrices carry also a time-index, i.e.
where the subscripts in tr c,d,t , det c,d,t indicate that these operators act on color, Dirac and time indices. C s 0 describes a closed spatial path of length s 0 (where, in contrast to (3.1), backtracking is now allowed) and the matrixM Cs 0 is given by the (ordered) product of the matrices S(1 − κ T ) −1 along that spatial path (see figure 16), x i+1 is understood to be a matrix with color, Dirac and time indices.
The final form of (3.78), which could be called a spatial loop expansion, is particularly useful if one is aiming towards including spatial fermion hoppings into the mean-field treatment. In this case one can just restrict the product over s 0 in (3.78) to the factors with s 0 = 0, 2 but still get the full interaction between nearest-neighboring sites. In contrast, in the original spatial hopping expansion (3.77), a nearest-neighbor truncation would limit the accuracy to O κ 2 s . In order to see how the presence of spatial hoppings changes the properties of the static-quark system, it is sufficient to take just the lowest order terms in κ s of (3.78) into account. These terms lead, after integrating out the spatial links, to corrections up to order O κ 2 s a nt f (β) (see appendix B). Using this improved effective nearest-neighbor Polyakov loop action within the mean-field framework described above in section 3.4, we can for example compute the average Polyakov loop for a mass degenerate two-flavor system. 3. and the one from section 3.4, which neglects fermion hopping completely (dotted red curve).
As can be seen, allowing for spatial fermion hopping in the effective theory shifts the value of µ where the average Polyakov loop is maximal closer to the corresponding value found in full LQCD. With increasing approximation order also the magnitude of the average Polyakov loop gets closer to the LQCD result. In the right-hand part of figure 17 we tried to compensate for the truncated spatial nearest-neighbor interaction in the effective theory by increasing the value of the inverse gauge coupling β. Setting β = 6.035 in the effective theory, the average Polyakov loop obtained by mean-field with the O κ 2 s a nt f (β) effective action from appendix B then almost coincides with the corresponding LQCD result for β = 5.0.
Finally, figure 19 shows the average sign (3.21) as a function of µ, again for a massdegenerate (κ = 5.0) two-flavor system of spatial size V = L 3 = 4 3 , computed by meanfield, using our improved effective action and However, the value of the hopping parameter κ = 0.15 used in figures 17 and 18 is already almost too large to be used in a calculation based on a O κ 2 s a nt f (β) truncated nearest-neighbor interaction. For larger values of κ or lower temperatures (i.e. larger N t ), the deviation of the mean-field result from the Monte Carlo data becomes larger and it is no longer possible to match the data by adjusting the inverse gauge coupling in the effective model as was done in the right-hand part of figure 17. This is not surprising as, concerning the physics, the O κ 2 s a nt f (β) truncation of the nearest-neighbor interaction term, det c,d,t 1 − κ 2 sM C 2 , implies that mesons and baryons (the low energy degrees of freedom of the theory) remain essentially static and cannot undergo spatial hoppings. But when lowering the temperature and/or the quark mass, it is precisely the hopping of these low energy degrees of freedom that becomes more important. From a more technical point of view, looking at the expansion of the fermionic nearest-neighbor interaction term,

JHEP02(2016)051
where the c k (A) for a rank-n matrix A are defined by the characteristic polynomial of A, it becomes clear that the κ 2 s term does not necessarily give the dominant contribution to (3.80) as the c k M C 2 can become rather big (especially if n t is large) and to guarantee that the κ 2 s -term dominates the expansion, one would have to choose κ s much smaller than we did.
To incorporate spatial meson and baryon hopping into the mean-field calculation, one should include at least the κ 6 s -terms in the expansion of (3.80). Although this is probably feasible, already at order κ 4 s the computation by hand of the corresponding terms in the effective action becomes a delicate issue and it would be desirable to automate the combinatorics required to get the correct dependency on a f (β).
In order to demonstrate that the deviation of the mean-field result from the Monte Carlo data in figures 17 and 18 is merely due to the truncation of the nearest-neighbor interactions and not caused by our mean-field method itself, we try to reproduce figure 6 of [6], which shows L and L * for a single flavor system with κ = 0.01, N t = 200 and β = 5.7, obtained once by a Monte Carlo and once by a Complex Langevin simulation of the effective model derived in [6]. Our mean-field action is obtained by applying (B.43) to the partition function (B.40), which is of order κ 2 s in κ s , while for [6, figure 6] also terms proportional to κ 4 s have been taken into account. However, as κ is rather small, we should nevertheless get a result comparable to [6, figure 6] if our mean-field method works correctly. In figure 20, left, we show [6, figure 6] superimposed with the corresponding result from our mean-field calculation, where also for this single flavor system, we determined the mean-field value L within the phase-quenched system and then used reweighting to obtain L and L * in the phase-unquenched case. As the figure shows, our mean-field results (dotted lines) match remarkably well the data obtained by full Monte Carlo and complex Langevin simulations. On the right hand side of figure 20, we also compare the mean-field result obtained with the O κ 2 s truncated effective fermion action with the corresponding result obtained with the "full" (O κ 2 s a Nt f (β) truncated) effective fermion action (B.33). As can be seen, even for the small value of κ = 0.01 and large N t = 200, the higher order terms in the effective fermion action lead to a clearly observable difference in the result.

Summary
After a short introduction to the sign problem, we have shown that the half-filling point, i.e. the value of the chemical potential where the fermion density assumes half its maximal value, has similar properties in LQCD as in e.g. the fermionic Hubbard model used in solid states physics: the sign problem almost disappears at half-filling and the system possesses an apparent particle-hole symmetry. We traced the origin of these properties by analyzing the heavy-dense fermion determinant and the corresponding strong coupling partition function. However, in contrast to the Hubbard model, where the lattice and the half-filling state are physical, in LQCD the lattice is just a regulator that has to be removed in order to extract continuum physics and the state corresponding to half-filling of that lattice is therefore of no physical relevance.
As the half-filling state in LQCD and the corresponding maximum in the average sign are just lattice artifacts, it seems sensible that the first appearance of a minimum in the average sign (as a function of increasing chemical potential) indicates the point where lattice artifacts start to become dominant and measurements of observables taken at larger values of the chemical potential should therefore be considered with caution.

B Higher order corrections to the effective Polyakov loop action
In this section we derive an effective nearest-neighbor Polyakov loop action that includes not just the nearest-neighbor effects coming from the gauge field action, but also those coming from the fermion determinant. We start with the single flavor partition function which, using the spatial loop expansion (3.78) for the fermion determinant and the character expansion (3.61) for the gauge field Boltzmann factor, can be written as: where x, y are nearest-neighboring sites spanning a path C 2 . Equation (B.3) would capture the full spectrum of nearest-neighbor interactions but is unfortunately still rather complicated. We simplify the expression further by expanding the determinants det c,d,t 1 − κ 2 sMx,y to order κ 2 s and by considering the gauge field only in the fundamental representation and only along temporal plaquettes. With these simplifications, the two products in (B.3) turn into the following expression: The Boltzmann factor for the desired effective action should then coincide with (B.4) up to order O κ 2 s a nt f (β) (the single site/non-interaction terms coming from the factor Det(1 − κT ) in (B.2) are not considered as part of the effective action).
To find such an action, we start by writing out tr c,d,t M x,y with y = x + i and i ∈ {1, 2, 3}: · tr c ( U i (x, x 4 ) P y (x 4 , y 4 ) 1 + (2 κ e µ ) nt P y (y 4 ) −1 Integrating over the spatial links in (B.6) at order O a 0 f (β) leads to the constraint y 4 = x 4 (see figure 21) as for SU(3) integrals, we have: According to the definition of P x (x 4 , y y ) given below eq. (3.29), we then have P x (x 4 , x 4 ) = 1 and find therefore: which is the same as eq. (2.28) of [6], as can be seen by noting that: and tr c 2 κ e −µ nt P † To obtain non-zero contributions from terms in (B.6) for which y 4 = x 4 , we have to take plaquette terms into account. To see in more detail how this works, let us as an example consider the top left diagram in figure 22. It shows a contributing diagram to the first term inside the large bracket in the sum over x 4 , y 4 in eq. (B.6) for the case where y 4 − x 4 = 2. If we would integrate out the spatial links in the corresponding term in (B.6) without taking into account plaquette terms, the result would just vanish due to the left-hand identity in (B.7). But allowing for plaquette terms coming from the non-trivial gauge-field part JHEP02(2016)051 in (B.4), we can form the two left-most diagrams on the second and third row of figure 22, and get for example for the one on the second row, a term like and Integrating over the three spatial links U i (x, x 4 ), U i (x, x 4 + 1) and U i (x, x 4 + 2) then yields where, as per definition (see section 3.2.3), and P x (x 4 + 2, x 4 ) = U 4 (y, x 4 + 2) . . . U 4 (y, n t − 1) U 4 (y, 0) U 4 (y, x 4 − 1) , (B.16) we can simplify U † 4 (y, x 4 + 1) U † 4 (y, x 4 ) P y (x 4 , x 4 + 2) = 1 , and to get either the identity or a closed Polyakov loop. With this, we finally obtain for eq. (B.14): In a similar way we can proceed with the other terms in (B.6) and find, after grouping them according to the power of their factors of a f (β): 2 mod(y4−x4,nt) tr c 1 + (2 κ e µ ) nt P y −1 tr c 1 + 2 κ e −µ nt P † + (a f (β)) mod(x4−y4,nt) tr c (2 κ e µ ) nt P y 1 + (2 κ e µ ) nt P y −1 tr c 1 + (2 κ e µ ) nt P x −1 Note that all the traces in (B.20) are time-independent. The time-dependency of the individual terms in the sum is only in the different exponents of a f (β), (2 κ) 2 and (2 κ) −2 which form essentially simple geometric sequences. The (finite) sum can therefore easily be carried out and yields: (B.21) For x 4 = y 4 , there is a third class of diagrams, which are always of order κ 2 s a nt f (β),

JHEP02(2016)051
depicted in the last row of figure 22. If we take as an example again the left-most diagram on that row, we obtain by integrating out the spatial links, using the identity a term of the following form: we can simplify: , if x 4 = y 4 , (B.25) and (B.23) becomes for x 4 > y 4 : − a nt f (β) 2 (2κ e µ ) nt tr c P y tr c P y 1 + (2κ e µ ) nt P y −1 − tr c P 2 y 1 + (2κ e µ ) nt P y and for x 4 < y 4 : − a nt f (β) 2 (2κ e µ ) nt tr c P † y tr c 1 + (2κ e µ ) nt P y −1 − tr c P † y 1 + (2κ e µ ) nt P y −1 tr c P x tr c P x 1 + (2κ e µ ) nt P x −1 − tr c P 2 x 1 + (2κ e µ ) nt P x (2κ e µ ) nt tr c P y tr c P y 1 + (2κ e µ ) nt P y −1 − tr c P 2 y 1 + (2κ e µ ) nt P y −1 tr c P x tr c P x 1 + (2κ e µ ) nt P x −1 − tr c P 2 x 1 + (2κ e µ ) nt P x
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.