Sign problem in finite density lattice QCD

The canonical approach, which was developed for solving the sign problem, may suffer from a new type of sign problem. In the canonical approach, the grand partition function is written as a fugacity expansion: $Z_G(\mu,T) = \sum_n Z_C(n,T) \xi^n$, where $\xi=\exp(\mu/T)$ is the fugacity, and $Z_C(n,T)$ are given as averages over a Monte Carlo update, $\langle z_n\rangle$. We show that the complex phase of $z_n$ is proportional to $n$ at each Monte Carlo step. Although $\langle z_n\rangle$ take real positive values, the values of $z_n$ fluctuate rapidly when $n$ is large, especially in the confinement phase, which gives a limit on $n$. We discuss possible remedies for this problem.

The canonical approach, which was developed for solving the sign problem, may suffer from a new type of sign problem. In the canonical approach, the grand partition function is written as a fugacity expansion: Z G (µ, T ) = n Z C (n, T )ξ n , where ξ = exp(µ/T ) is the fugacity, and Z C (n, T ) are given as averages over a Monte Carlo update, z n . We show that the complex phase of z n is proportional to n at each Monte Carlo step. Although z n take real positive values, the values of z n fluctuate rapidly when n is large, especially in the confinement phase, which gives a limit on n. We discuss possible remedies for this problem. 1. Sign Problem and the Canonical Approach A lattice QCD simulation is a first-principles calculation, and this makes it possible to study the quark/hadron world using a nonperturbative approach. The basic formula is the Feynman path integral form of the grand partition function: where µ is the chemical potential, T is the temperature,Ĥ is the Hamiltonian,N is the number operator, det ∆(µ) is the fermion determinant, and S G is the gluon kinetic energy.
In this paper, we consider the two-flavor case: N f = 2.
To explore the finite density QCD, we consider finite µ regions. However, when µ takes a nonzero real value, the fermion determinant becomes a complex number. This is problematic, because in the Monte Carlo simulations, we generate the gluon fields with the probability and if the fermion determinant is complex, we are in trouble. In principle, we may write det ∆ = | det ∆| exp(iφ), perform the Monte Carlo update with | det ∆|, and push the phase exp(iφ) into an observable. However, the fermion determinant may have the form det ∆ = exp(−f V /T ), and the phase fluctuation Imf V /T becomes large as the volume becomes large, so this does not work in practice.
Recently, the canonical approach, or the fugacity expansion, has attracted much attention as a candidate for solving the sign problem [1][2][3][4][5][6][7][8][9][10]. In the canonical approach, the grand canonical partition function is expressed as a fugacity expansion: where ξ = exp(µ/T ) is the fugacity. Both Z G and Z c are functions of the volume,V , which we abbreviate in the arguments. The inverse transformation is In Eq. (3), the canonical partition functions, Z C , do not depend on µ, and Eq. (3) works for real, imaginary, and even complex µ. When the chemical potential is pure imaginary, µ = iµ I , the fermion determinant is real, and in those regions, we can construct Z C from Z G . After determining Z C in this way, we can study real physical µ regions using formula (3).
2. Calculation of Z C (n) In order to obtain the canonical partition functions, Z C (n, T ), first we expand the fermion determinant: Then, Here, · · · 0 is an expectation value at µ 0 . One can assign any pure imaginary value to µ 0 . There are various ways to obtain z n : (1) Direct evaluation of det ∆ [4].
(3) Winding number expansion [12,13] In this letter, we employ method (3), in which the fermion determinant det ∆ is expanded as a series of the hopping parameter κ, and the diagrams are classified and packed with respect to the fugacity power: An algorithm that describes how to obtain Eq. (8) from Eq. (7) is given as "Algorithm 1: Winding Numbers via Hopping Parameter Expansion" in Ref. [12]. In Fig. 1, we show some of the winding diagrams of Eq. (7), which contribute W 0 , W 1 , and W −1 . In Fig. 2, the magnitude of the winding numbers |W k | is shown as a function of k, for M max = 300, 600, 840, and 1120 in the confinement temperature region (T /T c = 0.93).

Result
The lattice QCD simulations reported here were performed at the Far Eastern Federal University on Vostok-1, which consists of 10 nodes (2 × Intel E5-2680 v2, 64 GB RAM; 2 × Nvidia Tesla K40X Kepler). Its LINPACK performance is 23.52 TFlops. In our code, the clover Dirac operator performance was 76.9 GFLOPS (53.4% of the peak on the GTX 980).

3/7
The lattice size N x N y N z × N t is 16 3 × 4, and β and the hopping parameter are chosen to ensure m π /m ρ = 0.80 [14].
In Fig. 3, we show the complex phase of z n for several configurations. The first observation is that they are approximately proportional to n. In Ref. [6], the authors presented the relationship between the winding number expansion and the canonical partition functions, as follows: Here, and we use the relation Then, we obtain the Fourier transform of Eq. (10) to get Z C in Eq. (4) 1 , Using the integral representation of the modified Bessel function, the lowest order of Eq. (13) reads This is proportional to z n in each configuration. In this lowest order, we see that z n has a phase coming from e inδ1 , which is proportional to n.
The phase of z n is approximately proportional to n. In the following, we parametrize the phase as nδ. In the deconfinement phase, the slope is small, namely, z n are nearly real, while in the confinement phase, the slope is large, sometimes crossing ±π.
The distribution of δ is shown in Fig. 4, and the very different behavior in the confinement and deconfinement phases is apparent. Above T c , the phase of z n is almost zero, while below T c , the phase fluctuates significantly. For example, if δ = 0.1, nδ > π for n > 31. In other words, for large n, the real part of z n fluctuates between positive and negative values in the confinement phase. After averaging many configurations, the average of z n should be real positive. But when n is large, we suffer from a "sign problem". In Fig. 5, we show the behavior of the fermion determinant in the regions in which the chemical potential is pure imaginary. The data are evaluated by the reweighting method at µ 0 I /T = 0, 2π/3, and 4π/3, in order to recover Roberge-Weiss symmetry.

Conclusion
In this letter, we reported that there is a hidden sign problem in the canonical approach, namely, z n having a complex phase. This was observed in Refs. [15] and [12]. We have confirmed it and found that it produces positive/negative cancellation in the confinement phase. We studied the distribution of the phase, δ, as shown in Fig. 4. This provides an approximate upper bound on |n| of Z n : |n| < π/ |δ| . The arguments from Eqs. (10) to (15) tell us that the phase of the winding number W n determines the phase of z n . Although each W n contributes z l , the largest effect comes from W ±1 . These consist of many diagrams, but the ones with the largest contribution are the Polyakov lines, L and L † . Therefore, the phase of the Polyakov lines produces the phase of z n . Values of L and L † scatter around the real axis, and at lower temperatures, |L| is small, which results in the large phase.
It is known that, in the grand canonical approach, the complex Polyakov line contributes the phase of the determinant. In the canonical approach, we can pinpoint the origin of the sign problem. In order to overcome the sign problem, it is important to pursue the following: (1) Study the lattice size and the quark mass dependence. In particular, we must investigate whether this new sign problem increases or decreases in severity as the lattice volume increases. (2) Find a solution to reduce this sign problem. The canonical approach can be combined with the Lefschetz thimbles method or other practical methods [16], [17].