Nonanalyticity, sign problem and Polyakov line in Z3-symmetric heavy quark model at low temperature: Phenomenological model analyses

The nonanalyticity and the sign problem in the Z3-symmetric heavy quark model at low temperature are studied phenomenologically. For the free heavy quarks, the nonanalyticity is analyzed in the relation to the zeros of the grand canonical partition function. The Z3-symmetric effective Polyakov-line model (EPLM) in strong coupling limit is also considered as an phenomenological model of Z3-symmetric QCD with large quark mass at low temperature. We examine how the Z3-symmetric EPLM approaches to the original one in the zero-temperature limit. The effects of the Z3-symmetry affect the structure of zeros of the microscopic probability density function at the nonanalytic point. The average value of the Polyakov line can detect the structure, while the other thermodynamic quantities are not sensible to the structure in the zero-temperature limit. The effect of the imaginary quark chemical potential is also discussed. The imaginary part of the quark number density is very sensitive to the symmetry structure at the nonanalytical point. For a particular value of the imaginary quark number chemical potential, large quark number may be induced in the vicinity of the nonanalytical point.


I. INTRODUCTION
Study of the quantum chromodynamics (QCD) phase structure at finite temperature T and quark chemical potential µ is one of the most important subjects in particle and nuclear physics, astrophysics and cosmology. Nowadays, the first-principle nonperturbative calculation, the lattice QCD (LQCD) simulation has been almost established at µ = 0. However, for µ = 0, LQCD has a famous sign problem and is very difficult to be done correctly. An effective action obtained after the integration over the quark fields is complex and the numerical simulation such as the Monte Carlo simulation is very difficult, since we can not construct a proper probability density function. Several methods were proposed so far to circumvent the sign problem; namely, the reweighting method [1], the Taylor expansion method [2,3] , the analytic continuation from imaginary µ to real µ [4][5][6][7][8][9], the complex Langevin simulation [10][11][12][13][14][15] , the Picard-Lefschetz thimble theory [16][17][18][19], and the path optimization method [20,21]. Particularly for the case of µ/T > 1 , our understanding of the QCD phase diagram is still far from perfection.
It was also suggested that the sign problem may be weaker in the Z 3 -symmetrized QCD than the original one [22]. Due to the effects of dynamical quark, the Z 3 symmetry which exists in the pure gluon theory and is related to the quark confinement is explicitly broken. However, in the symmetric three flavor QCD, the Z 3 symmetry can be restored by introducing imaginary isospin chemical potential with the absolute value i 2 3 πT . In this paper, we call the Z 3 -symmetric QCD "Z 3 -QCD" [22][23][24][25][26][27][28][29][30][31][32]. In Z 3 -QCD and its effective models, the sign * kounoh@cc.saga-u.ac.jp † Kashiwa@fit.ac.jp ‡ hirakida@izumi.ac.jp problem is expected to be weaker than the original ones, since the number of configurations of which the effective action are real increases by symmetrizing the theories. In fact, the Z 3symmetric three-dimensional three state Potts model has no sign problem [33]. In the Z 3 -symmetric effective Polyakovline model (EPLM), the sign problem remains, but it is much weaker than in EPLM without Z 3 -symmetry [34]. (In this paper, we call the Z 3 -symmetric EPLM "Z 3 -EPLM". ) Figure 1 shows the schematic phase diagram obtained by using Z 3 -EPLM with the reweighting method [34]. This diagram is also consistent with the one obtained by using the Z 3 -symmetric three-dimensional three-state Potts model with no sign problem [33]. In Z 3 -EPLM, the sign problem happens only nearby the line µ = M at low temperature. Since EPLM is the model of heavy quarks, the chiral symmetry restoration is not expected. Hence, it is naturally expected that this sign problem is simply related to the Fermi sphere formation at µ = M . In this meaning, here we call this sign problem "a trivial sign problem". To detect anomalous phenomena, we need to remove or weaken it.
It is well known that the imaginary chemical potential is transformed into the change of the temporal boundary condition of quark fields by redefining the quark fields [35]. The temporal boundary condition is irrelevant in the zero temperature limit, namely β = 1/T → ∞. Hence, Z 3 -QCD (and its effective model) approaches to the original QCD (and effective model). In fact, in Ref. [22], it was shown that the phase diagram of the Z 3 -symmetric Polyakov-loop extended Nambu-Jona-Lasinio (PNJL) model coincides with the one in the original PNJL model [36][37][38][39][40] in the zero-temperature limit. However, this limit may be nontrivial at finite µ since nonanalyticity occurs at zero temperature due to the Fermion sphere formation and this phenomena itself is also related to the change of the boundary condition and the zeros of the grand canonical partition function. Furthermore, in Z 3symmetric theory, the expectation value of the Polyakov line (loop) vanishes due to the exact Z 3 -symmetry, while it can be finite in the original model. It is nontrivial whether the Polyakov line coincides or not in two models in the zerotemperature limit. (Note that the expectation values of the absolute value of the spatial average of the Polyakov line can be finite in Z 3 -symmetric model and is used to analyze the confinement-deconfinement transition [33,34]. ) In this paper, using the heavy quark model, we study phenomenologically how the Z 3 -symmetric model approaches to the original one and weaken the sign problem. We also examine how the Z 3 -symmetry affects the nonanalyticity at zero temperature. This paper is organized as follows. In Sec. II, we study the analyticity in Z 3 -symmetric heavy quark model in the relation to the zeros of the grand canonical partition function Z [41,42] and the boundary condition of the quark field. The zeros structure of Z of the free lighter fermion gas with spatial momentum is also discussed. In Sec. III, using Z 3 -EPLM in the strong coupling limit, we examine how Z 3 -EPLM approaches to the original EPLM and weaken the sign problem in the zero temperature limit. Relation among the Z 3 symmetry, the zeros of the probability density function and the sign problem is discussed. It is shown that the Polyakov line at the nonanalytic point can detect the symmetry structure of the zeros, while the other quantities are not sensible to the structure. The effect of the imaginary quark chemical potential and the Roberge-Weiss (RW) periodicity [35] at low temperature limit is also discussed. Section IV is devoted to a summary.

II. NONANALYTICITY IN FREE HEAVY QUARK MODEL AT ZERO TEMPERATURE
A. EOS of Free fermion at zero-temperature In this subsection, we briefly summarize the nonanalyticity of the equation of states (EOS) of the free fermion at T = 0. For the free fermion with finite mass M , thermodynamical quantities vanish at T = 0, when µ < M . When µ ≥ M , the pressure P of the free fermion gas at T = 0 is given by where g is the fermion degree of freedom including the number of spin states. The fermion number density and its derivative with respect to µ are given by and Note that these quantities are continuous at µ = M . However, the third derivative of the pressure is divergent at µ → M + 0. Hence the pressure is nonanalytic at µ = M .

B. Free heavy quark model
In this subsection, we consider the free heavy quark model (FHQM) on the lattice with N s , N f = 3 and N c = 3 where N s , N f and N c are the number of the spatial sites, the number of flavor and the number of color , respectively. Quarks in the heavy mass limit have no spatial momentum and the energy of them is always equal to their mass M f (f = u, d, s). Hence, the grand canonical partition function is given by where N = 2N c N s , and M f and µ f are the mass and the chemical potential for the f quark. In (5), the antiquark contributions are included since they are important for the reality of Z when imaginary chemical potential is introduced, although the contributions vanish in the limit M f , µ f → ∞. If we put µ f = M f + i(2k + 1)πT with an integer k, we obtain Z = 0.
When T approaches to zero, the location of the zeros of Z approaches a real value µ f = M f . Hence, in the analogy of the famous Lee-Yang theorem [41,42], the (dimensionless) pressure is expected to be non-analytic at µ = M f when β → ∞. (For the application of the Lee-Yang theorem to the QCD phase transitions, e.g., see [43] and references therein. ) Note that, beside the infinite limit of the spatial volume N s , here we also take the infinite limit of the imaginary time (τ ) length β. It is known that, by the redefinition of the quark filed the imaginary chemical potential µ I = iθT can be transformed into the twisted temporal boundary condition Hence, if θ = (2k + 1)π, the quark boundary condition is a periodic boundary condition. It should be noted that, except for the singular point µ = M f + i(2k + 1)π, the twisted model approaches to the original model in the limit β → ∞, since the boundary condition becomes irrelevant in this limit.
(Since θ has a trivial periodicity 2π, in this paper, we restrict θ in the region (−π, π] for simplicity. ) Hereafter, we consider the flavor symmetric case, namely, M u = M d = M s = M , unless otherwise mentioned. The derivative of P with respect to µ, namely, the (dimensionless) number density is given by and it is clear that n q is nonanalytic at µ = M and T = 0. However, in the actual numerical simulations such as EPLM, we can not put β = ∞. Hence, n q and its derivatives are continuous functions of µ as seen in Figs. 2 and 3. Near the point µ = M , the quark number density n q increases monotonically. The second derivative n ′′ q = ∂ 2 nq ∂µ 2 has maximum and minimum around the point µ = M . When M/T becomes larger, n q increases more rapidly and the absolute values of the maximum and minimum of n ′′ becomes larger, while the width of the peaks becomes narrower.
Figs. 4 and 5 are the same as Figs. 2 and 3, respectively, but for θ = π. As Re(µ) approaches to M , n q and n ′′ q diverge. The region of the divergent behavior be narrower as M/T be larger. Hence, it is expected that the results with θ = π approach to those with θ = 0 except for the point of Re(µ) = M . Here, we only show the results of the odd derivatives of the pressure P with respect to the chemical potential µ. As is seen in the next section, in EPLM, these odd derivatives are related to the sign problem around µ = M .
As is seen in the previous subsection, in the case of the free fermion with smaller mass and spatial momentum, the number density itself is continuous at µ = M . It seems that the effects of the spatial momentum make the transition smoother. At    finite temperature, the pressure of the free fermion gas is given by Then, the grand canonical partition function is given by is zero at µ = (n∆p) 2 + M 2 + iπ. Zeros of (1 + ) depend on the absolute value p = n∆p of the spatial momentum of fermions. The set of zeros of Z forms a continuous structure. This structure of zeros of Z may smoothen the transition. In Fig. 2, the quark number density n q in FHQM with M u = 0.5M , M d = M and M s = 1.5M is shown for M/T = 10. We see that n q increases slowly as µ increases. For the infinite flavor quarks with mass (n∆p) 2 + M 2 with infinitesimal degree of freedom, gn 2 ∆p 3 4π 2 , it can be expected that n q increases slowly, even if the limit M/T → ∞ is taken. In this meaning, the result in FHQM with nonsymmetric flavors mimics the free quark model with lighter mass and the spatial momentum.

C. Z3-symmetrization
In the symmetric three flavor quark model, we consider the case where In this paper, we call this setting "Z 3 -symmetrization" and Z 3 -symmetric FHQM "Z 3 -FHQM". Since the additional imaginary chemical potential iθ f T is the isospin chemical potential rather than the quark chemical potential, it can not be included in the definition of µ. It should be noted that Z is real at µ f = µ + iθ f T for real µ and is not zero The setting (12) of chemical potential is related to the socalled Z 3 symmetry. In fact, by the redefinition of the quark field q f , the imaginary part of the chemical potential µ f can be transformed into the temporal (τ ) boundary condition (13) where e −iθ f is an element of the Z 3 group. In QCD, the Z 3transformation changes the boundary condition of quark field by the factor of the Z 3 group element. Hence, the Z 3 symmetry which exists in the pure gluon theory and related to the quark confinement is explicitly broken. However, in the symmetric three flavor QCD, the Z 3 symmetry is restored by use of the Z 3 -symmetrization (12) [24]. It should be noted that the Z 3 -symmetric theory is expected to approach to the original one in the limit β → ∞, since the boundary condition is not relevant in the limit.
In Z 3 -FHQM, the partition function Z becomes zero at following three points, Hence there are three zeros of Z in the complex µ plane. However, these zeros correspond to the same nonanalyticity at µ = M when the zero temperature limit is taken. Therefore, Z 3 -symmetric FHQM has the same information of the nonanalyticity of the original FHQM at zero temperature, although the two models are different each other at finite temperature. It is known that the Z 3 -symmetrization enhances the confinement-like structure and a quark behaves as the particle with the mass 3M rather than that with M . Hence, the effects of the Z 3 -symmetrization resembles the ones of the increases of M/T . In fact, the partition function of Z 3 -FHQM is given by Hence, Z 3 -FHQM has the same properties as the ordinary FHQM with one flavor, threefold mass and threefold chemical potential. If we define the baryonic chemical potential µ B = 3µ, the zero of Z locates only at µ B = 3M + iπ. The quark number density is given by The factor 3 complements the flavor decreasing from 3 to 1 in Eq. (15). Hence, the number density in Z 3 -FHQM coincides the one in the ordinary FHQM with three flavor, threefold mass and threefold chemical potential. The similar correspondence is seen in n ′′ M 2 = ∂ 2 nq ∂(µ/M) 2 . In Figs. 2∼5, the Re(µ)-dependence of n q and n ′′ q of the Z 3 -FHQM are shown for θ = 0 and π. In all case, the perfect coincidence between Z 3 -FHQM with M/T = 10 and FHQM with M/T = 30 is seen. When M/T is fixed and θ = 0, n q increases more rapidly and the absolute values of the maximum and minimum of n ′′ becomes larger in Z 3 FHQM than in FHQM, while the width of finite n ′′ becomes narrower. The localization of the peaks of the odd derivatives makes the sign problem weaker than that in EPLM.

III. EFFECTIVE POLYAKOV-LINE MODEL AT ZERO TEMPERATURE
A. Effective Polyakov-line model The grand canonical partition function of EPLM in temporal gauge is given by [11,15,34] where U x is the Polyakov line (loop) holonomy, the symbol i is an unit vector for i-th direction and L F is the fermionic Lagrangian density the concrete form of which will be shown later. The site index x runs over a 3-dimensional lattice. Large (small) κ in EPLM corresponds to high (low) temperature [44] in QCD. Although the relation between κ and T is not simple, we regard the case with κ → 0 and M/T → ∞ as the zerotemperature limit in this paper. The Polyakov line holonomy U x is parameterized as [15] with the condition ϕ r,x + ϕ g,x + ϕ b,x = 0. Instead of U x and U † x , here the phase variables ϕ r,x and ϕ g,x are treated as dynamical variables. The Haar measure part L H is given by [15] The (traced) Polyakov line (loop) P x and its conjugate P * x are defined as For the fermionic Lagrangian density with the flavor dependent quark chemical potential µ f and temperature T , we consider a logarithmic one of Ref. [15,34]: where M f is the quark mass. In the limit The reality of L F at µ = M f is related to the particle-hole symmetry [34,45]. In the zero-temperature limit, the antiparticle part of (26) vanishes. It should be also noted that, at µ = M f , there is no dependence on M f /T in EPLM with symmetric flavors, if the antiquark contributions can be neglected. Hence, the point µ = M is the fixed point where physical quantities are not changed when M f /T varies. Breaking of the flavor symmetry also breaks this invariance.
The (dimensionless) quark-number density n q is obtained by where N s is the number of the lattice spatial sites. As in the case of FHQM, Z 3 -symmetrization (12) can be done for EPLM with three symmetric flavors. In this paper, we call the Z 3 -symmetric EPLM "Z 3 -EPLM". The fermionic Lagrangian density of Z 3 EPLM is given by In the case of EPLM, the internal dynamical variables ϕ c,x are also threefold. Hence, this breaks the equivalence between Z 3 -EPLM and the ordinary EPLM with threefold mass and threefold chemical potential.
B. EPLM at κ = 0 For κ = 0, the partition function becomes simple, since the S G term vanishes. For large N s in which the periodic boundary condition is negligible, the integration over ϕ r,x and ϕ g,x can be performed independently for each site x. The partition function turns out to be where L = L H + L F , N s is the number of the spatial lattice sites, and z is the local partition function at one lattice site. It is known that the integral such as (29) can be evaluated analytically [45]. Hence, we call the equation "analytical" in this paper. However, the integral form is useful here, since we are interested in the mechanism of the sign problem. When L is not real, instead of L, we may use an approximate real Lagrangian L ′ for constructing the probability density function. Then, the approximate partition function reads (30) and the reweighting factor is W = Z/Z ′ = (z/z ′ ) Ns . When we put L ′ = Re(L), we obtain the reweighting (phase) factor in the phase quenched (PQ) approximation. In this paper, we call the reweighting method with PQ approximation "PQRW". We examine how PQRW works in EPLM. For the brief review of PQRW, see appendix A.
Using Eq. (29), the pressure P , the quark number density n q , the scalar density n s , the averaged value of the Polyakov line P x and its conjugate P x are given by respectively. These physical quantities are independent of N s , although Z depends on N s . The analytical form of these quantities are modified when we consider the boundary condition.
See Appendix B for detail. (It is easily seen that the modified results coincides with the equations above in the thermodynamical limit N s → ∞. ) In Ref. [34], it was found that the phase factor and physical quantities at small κ is very close to the ones at κ = 0. Hence, we use Eqs. (29)∼(35) as a phenomenological model for QCD with heavy quarks at low temperature. By use of these equations, we can discuss the sign problem from the results which are free from the sign problem. In the numerical calculations, we put M u = M d = M s = M unless otherwise mentioned. Figure 6 shows the µ-dependence of the phase factor W in PQRW. Since W depends on N s , we set N s = 10 3 , 20 3 and 30 3 . Around µ = M , W is small and the sign problem is serious in that region. However, due to the particle-hole symmetry, S F is real and W = 1 is always hold at µ = M . The small W indicates that the sign problem is serious when PQRW is used in simulations. When N s increases, the sign problem be more serious. However, the N s -dependence becomes small when N s is large. Hence, we set N s = 30 3 hereafter. Figures 7 shows the similar as Fig. 6 but the one with fixed N s . Roughly speaking, the sign problem is serious when |µ − M | < 5T . Hence, when M/T increases, the sign problem be weaker for fixed µ/M . (However, the situation may be different when we compare them for fixed µ. The phase factor W can be larger for lighter quark than for heavier one when we compare them at fixed µ. ) Figures 8 shows the result in EPLM without three flavor symmetry. When the strange quark mass M s is larger than the light quark mass M l , the sign problem becomes weaker, but the change is not so large. This indicates that the lighter quark dominates the sign problem. In Fig. 8, the results in Z 3 -EPLM are also shown. It can be seen that the Z 3 -symmetrization makes the sign problem weak drastically. In the case of Z 3 -EPLM with M/T = 100, the sign problem almost vanishes except for the vicinity of µ = M . In this case, from W itself, we see that nonanalyticity certainly happens just at µ = M .
Comparing Fig. 8 with Fig. 7, we see that the sign problem is somewhat weaker in Z 3 -EPLM with M/T = 10 than EPLM with M/T = 30. This is because, different from the number density, the partition function in Z 3 -EPLM is close to the one in EPLM with threefold mass and three fold chemical potential but one flavor. As well as the decrease of N s , the decreases of N f make the sign problem weaker. Figures 9 shows the µ-dependence of the quark number density n q . It is seen that n q abruptly increases at µ = M , when M/T is large. The results are antisymmetric with respect the point µ = M . This property is the consequence of the particle-hole symmetry. Due to the effect of the gauge field ϕ c,x , the equivalence between Z 3 -EPLM and EPLM with threefold mass and threefold chemical potential is slightly broken, but the difference be smaller when M/T becomes larger.
In the heavy quark model, the scalar density n s is almost the same as n q , since the effects of the spatial momentum and the vacuum fluctuations are absent and the antiquark contribution is negligible. If it couples to the quark field, it can make the quark mass smaller. In Fig. 10, the results in EPLM and Z 3 -EPLM in which quark mass changes from large mass M (= 10T ) to small one m(= T ) at µ/M = 1.1. Remember that W can be larger for lighter quark than the one for heavier quark when we compare them at fixed µ. Due to the change of the quark mass, the symmetry with respect to the line µ/M = 1 is broken. In EPLM and Z 3 -EPLM, the breaking of symmetry around µ = M may indicate a nontrivial change of the system. However, in QCD, such symmetry is not expected at the beginning. Hence, we should control the trivial sign problem caused by the formation of Fermi sphere anyway. Figures 11 shows the µ-dependence of the averaged value P x and P * x in EPLM. It is seen that the both quantities somewhat large in the vicinity of µ = M . P x has a maximum in the region µ > M , while P * x does in the region µ < M . This property is also a consequence of the particlehole symmetry. In Z 3 -EPLM, the expectation value of the Polyakov-loop vanishes due to the exact Z 3 -symmetry.
When all ϕ c,x vanish, P x and P * x becomes 1. The case corresponds to the ordered phase and the sign problem does not happen since all ϕ c.x vanish. However, the absolute values of the Polyakov line is far from 1 and the ϕ c,x fluctuates almost randomly. This causes the serious sign problem around µ = M . In Fig. 11, when M/T increases, P x and P * x becomes smaller at µ = M , but they do not change at µ = M since the point µ = M is the fixed point. Hence, in the limit M/T → ∞, the Polyakov line (and its conjugate) at µ = M and all the other quantities at any µ approach to the ones in Z 3 -EPLM, however, the Polyakov line (and its conjugate) at µ = M has different value in two models. It seems that the Polyakov line on the nonanalytical point can detect the difference of the boundary condition even in the zero-temperature limit.
It should be remarked that the existence of the fixed point at µ = M has an important role in the anomalous phenomena mentioned above. When the flavor symmetry is broken, the exact fixed point disappears. Figure 12 shows the µdependence of P x in EPLM with nonsymmetric flavors. In this case, the absolute values of these quantities are smaller than those in the symmetric flavor case. Furthermore, since µ = M f (f = u, d, s) is an approximate fixed point but not an exact one, the maximum values of these quantities decrease as M f /T increases. Hence, if we take the effect of the spatial momentum of quarks into account and take the zero temperature limit, the expectation value of the Polyakov line may vanish. (Althouh we do not show the result, the µ-dependence of P * x shows the similar tendency. )

C. Relation between sign problem and nonanalyticity at zero temperature
The local partition function z can be written as where the spatial index x is ommited for simplicity of the notation. The factor f = e −LF is given by where µ f,c = µ+iϕ c T in EPLM and µ f,c = µ+i(θ f +iϕ c )T in Z 3 -EPLM. This factor can be used as a microscopic probability density function in numerical simulation if the sign problem is absent. Figure 13 shows f /(|f | + ǫ) at µ = M in EPLM with M/T = 100, where ǫ is a positive infinitesimal constant. Note that, due to the particle-hole symmetry, f is real and is nonnegative at µ = M . Hence f /(|f | + ǫ) is 1 unless f = 0. The set of zeros forms a line structure. If one of Im(µ f,c ) is equal to (2k + 1)π, where k is an integer, f becomes zero at µ = M . This condition corresponds to the horizontal and vertical black lines at the edges in the figure. Furthermore, f is also zero when ϕ g = (2k + 1)π − ϕ r is satisfied. This condition corresponds to the black oblique lines in the figure.  Figure 14 shows the same as Fig. 13 but Z 3 -EPLM with M/T = 33.3 is used. Due to the Z 3 -symmetry, the zerostructure in the region where ϕ g = 2k−1 3 π ∼ 2k+1 3 π and 2l−1 3 π ∼ 2l+1 3 π with k, l = −1, 0, 1 is similar to the one of the EPLM result in the region where ϕ r = −π ∼ π and ϕ g = −π ∼ π.
The structure of zeros of f in Fig. 14 is Z 3 -symmetric but is not in Fig. 13. Hence, it can be said that the Polyako-line It should be noted that the zeros of f at µ = M themselves do not induce the sign problem. However, in the vicinity of µ = M , the situation is changed drastically. Suppose f = 0 at ϕ r = ϕ r0 and ϕ g0 when µ = M . Then, the absolute value of Im[L F (ϕ r,0 , ϕ g,0 )] may be still small at µ = M + ∆µ when |∆µ| is small enough. However, expanding L F (ϕ r , ϕ g ) at µ = M + ∆µ in term of ∆ϕ r = ϕ r − ϕ r,0 and ∆ϕ g = ϕ g − ϕ g,0 , we obtain L F (ϕ r , ϕ g ) = L(ϕ r,0 , ϕ g,0 ) +iL ′ F (ϕ r,0 ϕ g,0 )(∆ϕ r + ∆ϕ g ) where ′ denotes the differentiation with respect to µ. Since structure of the fermionic Lagrangian L F in EPLM is similar to the pressure of FHQM discussed in Sec. II, the odd coefficients in the expansion (38) have divergent behavior and Im(L F ) can be large at µ = M + ∆µ. This makes the sign problem serious in the vicinity of µ = M . Figure 15 shows There are area-like regions where Re(f ) is negative. These regions make |z| smaller than the quenched one z ′ and induce a serious sign problem. (Note that W will be almost zero in the large N s limit, even if |z| is slightly smaller than z ′ . ) Figure 16 shows the same as Fig. 15 but Z 3 -EPLM with M/T = 33.3 is used. As is the same as in Fig. 14, Z 3symmetric structure is seen also in this figure. The minimum value of Re(f )/(|f |+ǫ) is much larger than that in Fig. 15 and is not negative. Hence, the sign problem is not so strong in this case. When we set M/T = 100 in Z 3 -EPLM, Re(f )/|f | = 1 is realized anywhere, almost perfectly, the sign problem almost vanishes at µ = 0.98M .

D. Effects of imaginary quark chemical potential
When the imaginary quark chemical potential iθT is introduced, the reality of the grand canonical partition function Z and physical quantities is not ensured in general. Figure 17 shows the complex chemical potential dependence of Re(n q ) in EPLM with M/T = 100. Except in the neighbor of Re(µ) = M , the results for finite θ coincide with that for θ = 0. This is because the imaginary chemical potential, which is equivalent to the change of the boundary condition, is irrelevant in the zero temperature limit. In the neighborhood of Re(µ) = M the result for θ = π vibrates violently. The maximum value of |Re(n q )| is much larger than the degree of freedom of quarks, namely, 2N f N c = 18. In this figure, we show the result at 0.001 intervals in the horizontal axis. The singular behavior for θ = π depends strongly on the interval we use. When the interval is smaller, the maximum value of |Re(n q )| is larger. For θ = π, large quark number density can be induced at the singular point. This phenomenon happens also at Re(µ) = −M . Figure 18 shows the complex chemical potential dependence of Im(n q ) in EPLM with M/T = 100. When θ = 0 or π, Im(n q ) = 0. For θ = 5 6 π, Im(n q ) is finite only in the vicinity of Re(µ) = M . This is also because the imagi- nary chemical potential which induces the imaginary part of Im(n q ) is irrelevant in the zero temperature limit. Figures 19 and 20 shows the µ-dependence of the number density in Z 3 -EPLM with M/T = 100. The number density n q has similar properties as the one in EPLM, but the singularity for θ = π is very sharp. We also observe the oscillating behavior in Im(n q ) for θ = 5 6 π. It is known that the introduction of the imaginary chemical potential rotates the Polyakov line in the complex plane. Hence, we use the modified Polyakov line Q x = e iθ P x instead of P x itself. Note that Q x has the Roberge-Weiss (RW) periodicity, namely, Q x (θ + 2 3 π) = Q x (θ) , but P x does not [46][47][48][49]. Figures 21 and 22 show the complex chemical potential dependence of Re( Q x ) and Im( Q x ) in EPLM with M/T = 100. In these figures, the same tendency is seen as in the case of n q . In Z 3 -EPLM, Q x is always zero due to the exact Z 3 -symmetry. Figures 23 and 24 show the θ-dependence of n q and Q x at µ = M , respectively. The RW periodicity is clearly seen in these figures. Note that the imaginary parts of n q and Q x are the indicators of the RW-transition [46][47][48][49]. Here the RW pe-  riodicity is smooth and this property is not changed by varying M/T since the point Re(µ) = M is the fixed point. Hence, it is expected that the RW transition does not occur even in the limit M/T → ∞. Of course, this may be natural since we have set κ = 0. If the interaction between gauge field is switched on, a nontrivial transition may happen. The study in EPLM with finite θ and nonvanishing κ at low temperature limit is an interesting problem in future. (The study on the Z N -spin model with the external complex field and the interaction between the spins can be found in Ref. [50]. In that case, the hard sign problem induced by the external complex field was found. ) Figure 25 shows the θ-dependence of n q at µ = M in Z 3 -EPLM. The RW periodicity is seen in these figures. Furthermore, in Im(n q ), the higher frequency mode with the period 2 9 π is clearly seen. This property is related to the Z 3symmetry. The θ-dependence of Im(n q ) is very sensitive to the Z 3 -symmetry structure at Re(µ) = M .

IV. SUMMARY
In this paper, we have studied the non-analyticity and the sign problem in the Z 3 -symmetric heavy quark model at low temperature and examined how the Z 3 -symmetrized models approach to the original ones in the zero temperature limit. For the free fermion quark model (FHQM), the nonanalyticity at µ = M is related to the existence of zeros of the grand canonical partition function Z at finite temperature and complex chemical potential. By Z 3 -symmetrization, the zeros are threefold, but the Z 3 symmetric FHQM (Z 3 -FHQM) is equivalent to the original one with threefold quark mass and threefold quark chemical potential. Therefore, Z 3 -FHQM naturally approaches to the original one in the zero temperature limit.
We also examined the three flavor effective Polyakov-line model (EPLM ) with κ = 0. In Z 3 -symmetric EPLM (Z 3 -EPLM), the sign problem is drastically weaken in the low temperature comparing with the original EPLM. The Z 3 -EPLM also approaches smoothly to the original EPLM except for the nonanalytical point µ = M in zero temperature. At µ = M ,  the expectation values of the Polyakov line( and its conjugate) has different values in two models due to the existence or nonexistence of Z 3 -symmetry. The Polyakov line can detect the symmetry structure of the zeros of the microcanonical probability density function, while the other quantities are insensitive to the structure. This property is not changed by varying M/T , since µ = M is the fixed point in the flavor symmetric EPLM. However, the effects of the flavor symmetry breaking and the spatial momentum of quarks may break this property and the expectation value of the Polyakov line may vanish even in the original EPLM when M/T → ∞.
The effects of the imaginary chemical potential iθT at low temperature was also studied. The physical quantities at finite θ coincides with those at θ = 0 except for the neighborhood of Re(µ) = M . Hence the imaginary parts of the physical quantities can be induced only in the neighborhood. In the neighborhood of Re(µ) = M , the real parts of the number density and the modified Polyakov line vibrate violently. Large quark number density can be induced at the singular point. The θdependence of the imaginary part of the physical quantities at the nonanalytical point is affected by the symmetry structure of the microscopic probability density function.
It seems that the Z 3 -symmetrized theory is equivalent to the original one with larger mass at least except just on the nonanalytical point, when M/T is large enough. The trivial sign problem is expected to be weak in the Z 3 -symmetrized theory. Hence, to explore the low temperature property, we may use the Z 3 -symmetrized theory with smaller M/T instead of the original one. However, in LQCD, it is known that there is a nontrivial hard problem on early onset of quark number density at zero temperature. (See, e.g., Refs. [51,52] and references therein. ) Since Z 3 -QCD is expected to approach to the original QCD in zero temperature limit, this problem may also happen in the Z 3 -QCD when M/T is very large. However, the problem may not occur just below T c in Fig. 1 and we may use the probability density function in that region as the approximate probability density function to analyze the low temperature physics in the original QCD. Hence, the research of the lattice Z 3 -QCD at the intermediate temperature may be important. Such research is now in progress.
One of the simple approaches to the sign problem is the reweighting method. In this method, one can calculate the expectation value O ′ of the quantity O with a approximate weighting function F ′ (U ) which is real and nonnegative.
where U is the dynamical variables such as ϕ c,x in EPLM. The true expectation value O is given by where Z and W = Z/Z ′ are the true grand canonical partition function and the reweighting factor. When W is very small, the true expectation value has large errors due to the division by W in (A2). In actual calculations, the phase-quenched function is often used. This reweighting method is PQRW. In PQRW, W is also called as "phase factor".
Appendix B: Analytical representation of physical quantities in EPLM at κ = 0 with periodic boundary condition In the three flavor EPLM at κ = 0 with periodic boundary condition, the grand canonical partition function is given by where D is the dimension of the spatial space and L 3 s (= N s ) is the number of the lattice spatial sites. Note that Similarly, the partition function for an approximate Lagrangian L ′ , the pressure, the quark number density, the scalar density, the averaged values of the Polyakov line and its conjugate at κ = 0 are given by . (B7) In this case, not only the phase factors, but also the other physical quantities depend on N s . However, it can be easily seen that the effects of the boundary conditions vanish and the N sdependences of these thermodynamical quantities also vanish in the limit of N s → ∞.