QCD phase transition at real chemical potential with canonical approach

We study the finite density phase transition in the lattice QCD at real chemical potential. We adopt a canonical approach and the canonical partition function is constructed for Nf = 2 QCD. After derivation of the canonical partition function we calculate observables like the pressure, the quark number density, its second cumulant and the chiral condensate as a function of the real chemical potential. We covered a wide range of temperature region starting from the confining low to the deconfining high temperature; 0.65Tc ≤ T ≤ 3.62Tc. We observe a possible signal of the deconfinement and the chiral restoration phase transition at real chemical potential below Tc starting from the confining phase. We give also the convergence range of the fugacity expansion.


Introduction
The lattice QCD simulation provides a first-principle approach for the quark hadron world, and is expected to provide reliable information for the QCD phase diagram study. However, due to the notorious sign problem, its validity is limited; the fermion determinant DetD W (µ) appearing in the grand canonical partition function as a part of the probability weight is in general complex at real chemical potential µ, and the Monte Carlo calculation do not work.
In order to overcome the difficulty, several methods were proposed for studying the QCD at finite real chemical potential regions. One novel method was a multi-parameter reweighting [1], which initiated many studies [2]. Detailed analysis of the determinant ratios in this formula [3] revealed that the valid range is limited due to the overlap problem. The imaginary chemical potential method was proposed in refs. [4,5] and has become a useful tool; at the imaginary chemical potential the fermion determinant is not complex, but real. In ref. [6], it is shown that the polynomial fits to extrapolate the imaginary chemical potential regions to the real ones have large uncertainty for µ/T > 1, where T is the temperature. See ref. [7] for a nice review of these approaches and the overlap problem, and ref. [8] for more references on the two methods.
The grand canonical ensemble is thus a difficult subject to treat in lattice QCD because of the sign problem. The canonical partition function is related to the grand canonical one through the fugacity expansion and is known to be free from the complex action problem. The history of the canonical approach in finite density lattice QCD may have started with a development of a reduction formula for the Dirac determinant [9], which gives naturally a fugacity expansion of the Dirac determinant. A several studies with the staggered fermion are done along this line [10][11][12][13] and have been taken over by the Wilson fermion [14][15][16][17][18].
One of the problem of the canonical approach is that the Dirac determinant needs to be evaluated accurately and frequently, which is the heaviest part of the computation. One of the solution is to adopt the Taylor expansion [24]. In this paper we perform the fugacity expansion by a method of the temporal winding number expansion [20,21], which is represented as a hopping parameter expansion in temporal direction. By using the hopping parameter expansion the Dirac determinant can be evaluated with low cost and we can visit a wide parameter space in β which corresponds to 0.65T c ≤ T ≤ 3.62T c . The canonical partition function is evaluated for N f = 2 QCD both in the deconfinement and the confinement temperature regions. After derivation of the canonical partition function we study the chemical potential dependence of observables like the pressure, the quark number density and the chiral condensate. We observe a possible signal of the confinementdeconfinement and chiral restoration phase transition at real chemical potential below the critical temperature T c . A preliminary result has been reported in ref. [25].
This paper is organized as follows. In section 2 we briefly review the canonical approach. Our formulation of the fugacity expansion in terms of the hopping parameter expansion is discussed in section 3. After mentioning our numerical setup in section 4 our main results are given in section 5 for canonical and grand canonical partition functions. Section 6 is devoted for the chiral condensate in the grand canonical ensemble. A conclusion is given in section 7.

Canonical partition function
It is the well known fact that the grand canonical ensemble and the canonical one are equivalent each other. This is shown by a simple equation to relate the grand canonical partition function Z G (µ, T, V ) and the canonical Z C (n, T, V ) where ξ is the fugacity and this is a so called fugacity expansion formula. The inverse of this expansion is given by using the Cauchy's integral theorem Here we would like to emphasize an implicit assumption is made for the inverse transformation (2.2) to work. That is the grand partition function Z G (ξ, T, V ) should have singularities only at ξ = 0, ∞ corresponding to trivial ones at µ/T = ±∞. A phase transition point ξ c is not a singularity of the partition function but it is rather a zero of Z G (ξ, T, V ) (Lee-Yang zeros [26,27]) and does not affect the assumption. Now it is free to change the contour to a unit circle ξ = e iθ and the contour integral turns out to be a Fourier transformation [11] Z C (n, T, V ) =

JHEP02(2016)054
A standard Monte Carlo simulation is possible for the lattice QCD since the chemical potential is set to pure imaginary µ/T = iθ and the Boltzmann weight is real positive for two flavors. Eqs. (2.1) and (2.3) tell us that we construct canonical partition functions from the grand partition function at pure imaginary chemical potential, which are free from the sign problem, and that we can obtain the grand partition function at real chemical potential regions. Note that this is mathematically rigorous formula, and theoretically can be used at any large values of µ, as long as the right-hand side of eq. (2.1) converges. 1 However a difficulty of the sign problem is preserved unfortunately since there remains a frequent cancellation in plus and minus sign in a numerical Fourier transformation especially for large particle numbers. In order to get the canonical partition function accurately for large quark number n we need a very fine sampling of Z G (e iθ , T, V ) in θ. This requires a heavy computational cost because we need to evaluate the Dirac operator determinant for the QCD grand partition function. In this paper we shall solve this problem by a direct expansion of the Dirac determinant in terms of the fugacity.
In ref. [19], the canonical approach was tested on 4 4 lattice, and the authors further pursued the canonical approach [20,21], where the winding number in the expansion is limited up to six, and the Fourier transformation was done analytically. Then large baryon number regions are out of the scope.

Fugacity expansion by means of hopping parameter expansion
We consider the lattice N f = 2 QCD grand partition function given in the path integral form where we adopt the improved Wilson Dirac operator The term (3.5) is a so called clover term with the Sheikholeslami-Wohlert coefficient c SW [28], and F µν is the lattice field strength constructed from the standard clover-shaped combination of gauge links. Both the chemical potential e ±µa and the hopping parameter κ appear in front of the temporal hopping term simultaneously. The fugacity expansion of the determinant shall be given by using the hopping parameter expansion. The hopping parameter expansion originated, to our knowledge, from refs. [29] and [30]; the authors estimated dynamical fermion effects in the lattice QCD simulations. In ref. [31], the chemical potential was introduced in a hopping parameter expansion formulation, and a careful treatment of the small hopping parameter and large chemical potential leads to the high dense QCD model [32]. The hopping parameter expansion was combined with one link models to study the feature of finite density QCD [33]. Since the low order hopping parameter expansion cannot provide reliable information to analyze near the QCD phase transition, an improved hopping parameter expansion was proposed in ref. [34], where the temporal hopping is treated exactly, and the spatial hopping expansion was executed, where the method was employed in the complex Langevin simulation. Instead of the determinant we consider the hopping parameter expansion of its logarithm A non-trivial chemical potential dependence appears when the quark hoppings make a loop in the temporal direction. If one of the term has a n times winding loop in the positive temporal direction, the chemical potential dependence is e µanNt where N t is a temporal lattice length and this is nothing but e nµ/T = ξ n (See figure 1). Counting the winding number in temporal direction for each term in equation (3.6), we get the fugacity expansion where the coefficient w n (κ) contains O(κ m ) contributions of hopping term with m ≥ |n|N t . We shall explain the second step (3.7) in more detail. Let us consider a contribution from m quark hoppings given in a matrix form Q m (x, α; y, β), where x, y are site and α, β JHEP02(2016)054 represent spin-color index. Q m is expanded in terms of temporal hopping where k is a number of hopping in temporal direction. Negative value corresponds to that in negative direction. A positive and negative contribution cancels with each other in a single path and a total contribution turns out to be k. We have X starting from an initial condition Taking the trace the fugacity expansion coefficient is given by where sum should start with m = Min(4, N s ) for n = 0, with N s a spatial size of the lattice. The number 4 corresponds the minimal closed loop, the plaquette. Here we used a fact that the quark hopping should make a loop for the trace to survive. An explicit algorithm is given in ref. [35].
Exponentiating the summation we have a determinant The last step is to extract a coefficient z n (κ) in the following relation where factor two is for N f = 2 flavors. Then we get a fugacity expansion of the determinant [20,21] by means of the hopping parameter expansion Although the last step (3.13) is purely a mathematical issue no closed form is known for z n as a function of w n [20,21]. We use a numerical Fourier transformation to extract the coefficient

JHEP02(2016)054
The numerical Fourier transformation can be executed safely by using the multiple precision calculation [36] without any sign problem. This is because the fugacity dependence of the determinant or the partition function is given analytically once we evaluate the fugacity expansion coefficient w n . The fugacity expansion is done for gauge configurations generated at µ 0 = 0 or purely imaginary value. This procedure corresponds to the standard reweighting method The quantity we get in the Monte Carlo simulation is The canonical partition function is given by Now we are ready to calculate physical quantities in our canonical approach: 1. Generate gauge configurations at an imaginary chemical potential µ 0 by the standard HMC (Hybrid Monte Carlo) method.
3. Perform the Fourier transformation eq.(3.15), and calculate z n , which are averaged over the generated configurations. We have a canonical partition function Z C .
4. Using these averaged Z C , construct the grand partition function in eq. (2.1) at any real µ.

Numerical setup
We adopt the Iwasaki gauge action and the improved Wilson fermion action with the clover term. The number of flavors is set to N f = 2 with degenerate masses. The APE stout smeared gauge link [37] is used for those in the fermion action including the clover term. The number of smearing is four and the parameter is set to ρ = 0.1. The clover coefficient is fixed to c SW = 1.1 [38]. We adopt 8 3 × 4 lattice. A wide range of β is covered from high temperature β = 2.1 to low temperature 0.9. It seems to be that both the confining and deconfining regions are well covered, which is inferred from a behavior  Table 1. β and κ for the numerical simulation. m π /m ρ is also given. seems to be very near to the critical temperature T c as is shown the right panel of figure 2 using a fluctuation of phase of the Polyakov loop. A rough estimate of the temperature is given in the third column by using the gradient flow method, where we assume that β = 1.7 corresponds to the critical temperature for N t = 4. The gradient flow is used to get ratio of the lattice spacings for two neighboring β at zero temperature.
The hopping parameter is selected in order that the hopping parameter expansion works well, for which m π /m ρ turns out to be 0.7 -0.9 as is given in table 1. Maximal order of the hopping parameter expansion is set to 480 so that max winding number in temporal direction is 120. Number of independent configuration is 100 -600 as is shown in the table.

Canonical and grand canonical partition function
The first numerical result we get is the canonical partition function Z C (n, T, V ). We plot log |Z C (n, T, V )/Z C (0, T, V )|/(V T 3 ) as a function of the quark number n in figure 3. Data is shown only for quark number n which is multiple of three. The residual data do not JHEP02(2016)054 contribute to the fugacity expansion (2.1) because of the Z 3 symmetry. The partition function decays very rapidly with n and its behavior changes drastically between β = 1.9 (red) and 1.5 (green). The data above β = 1.9 is fitted very well by a quadratic function. But this is not the case below β = 1.5. This may correspond to a phase transition.
The plot range is fixed by using the d'Alembert's convergence condition which gives the convergence radius for the fugacity ξ. 2 The data are plotted in figure 4 for 1/3 log |Z C (n + 3)/Z C (n)|. If the thermodynamical system exist for QCD the Helmholtz free energy has a convexity in n, from which we can show a monotonic decrease of the quantity |Z C (n + 3)/Z C (n)|. We cut our data at n max where a monotonic decrease stops 3 indicated by vertical dotted lines in the figure. The horizontal dotted lines show maximal values of the fugacity expected to be within the convergence radius. Since the logarithm is taken, the line gives our applicable limit for the quark chemical potential µ/T . For example we can discuss physics safely at −3.2 < µ/T < 3.2 for β = 1.1 and −1.5 < µ/T < 1.5 for β = 2.1 with our method. We notice that n max is improved only by three even if we increase the maximal winding number twice as 120 → 240. This may be due to a fact that the convergence radius of the Taylor expansion (3.6) is unity. The second physical quantity is the grand partition function. By taking summation for −n max ≤ n ≤ n max in equation (2.1) with n multiple of three, we get the grand partition function for the real chemical potential. Here we notice that the canonical partition function Z C (n, T, V ) is a real positive quantity if a Hermitian transfer matrix exists for the lattice QCD. For example let us consider the Hamiltonian formalism of QCD in the continuum. The canonical partition function is given by a summation over every energy state with 2 For negative quark number we adopt |ZC (n − 3)/ZC (n)|. 3 This is due to a failure of the hopping parameter expansion at large quark number since a breakdown of the monotonic decrease is not seen when the exact formula is adopted for the fugacity expansion [15]. for a definition of the grand partition function. 4 Here we notice integers divisible by three contribute to the summation because of the Z 3 symmetry. From the partition function we get the pressure in the grand canonical ensemble
We plot the pressure normalized at µ = 0 in figure 5 as a function of the quark chemical potential. We observe that the pressure is very small at small chemical potential for low temperature below T c as is shown in the right panels of figure 5.  The third physical application is the k-th moment of the quark number operator

JHEP02(2016)054
In figure 6, the quark number density N /(V T 3 ) is plotted as a function of the real quark chemical potential µ/T . The data below T c is consistent with the hadron resonance gas prediction [39] N as is indicated by the solid line. This property is also understandable as a leading order JHEP02(2016)054 contribution from the fugacity expansion (5.5) The solid line in figure 6 is a plot of this contribution. The data is well described with (5.7) for T < T c since Z C (n, T, V ) decays rapidly with n as is shown in figure 3. The sinh function is not a good description for T ≥ T c . The number density rises linearly with µ/T near the origin as is indicated by the free quark gas model [40]. The second cumulant of the quark number density may be a good indicator of a phase transition. The results are plotted in figure 7. We observe peaks around µ/T ∼ 1 for every β. We also observe secondary peaks for T /T c = 0.72 (dark-green) and 0.83 (green) at larger µ/T , while only a single peak is observed at high temperature T ≥ T c . Although the secondary peaks are not shown for T /T c = 0.65 and 0.69 they exist at µ/T > 3. These peaks show a distinctive behavior against the quark number cut-off n max . The first peaks below T c are stable even if we vary n max by 20%. It would be conjectured that the first peaks in the low temperature side are the physical indication of the confinement-deconfinement phase transition. Here we notice that the quark mass is not the same for each β. Unfortunately we cannot exclude a possibility that these first peak corresponds to an onset phenomenon. We need a further study at lighter quark mass to have a firm conclusion. On the other hand, the second peaks at T /T c = 0.72, 0.83 and the first peaks above T c change their height and position with n max changes. These peaks would be artifacts due to the cut-off in quark number in the fugacity expansion.
The position of the first peak below T c is consistent with a physical expectation that the position moves to larger µ/T c as temperature decreases. We plot a rough estimate of the peak position in µ − T plane in figure 8. We set µ c = 0 at T = T c since β = 1.7 is very near to the phase transition point at µ = 0. We notice this is a preliminary plot since the quark mass is not the same between different β and no error estimation is performed. We observe two points at µ/T c ∼ 0.73. This is mainly due to an anomalous behavior of β = 0.9 data for which the peak position is rather low. This may be due to a large lattice artifact in such a strong coupling region.
The thermodynamic quantity as a function of the quark number is one of the most popular quantity in the canonical approach [13,17,21,24]. We notice the quantity plotted in figure 4 can be interpreted as an expression of the derivative in terms of the differential In refs. [13] and [17], an S-shape behavior was observed for N f = 4 and 3. In ref. [21], the authors see the S-shape behavior for N f = 4, but not for N f = 2 in the regions they searched, although there is some S-shape like structure.
Our results are plotted in figure 9. The data are very consistent with those of ref. [21], i.e., we see some structure, but cannot confirm the typical S-shape behavior. Signal is very weak and S-shape curve is obscure. We shall need more statistics to make the signal clear.

Hadronic observables
In our procedure, the fugacity expansion is based on the hopping parameter expansion. It may be possible to expand any hadronic operators in terms of the fugacity. We consider a JHEP02(2016)054 where O quark is a Wick contraction of a hadronic operator O in quark fields. For example, we consider the chiral condensate. It is easy to expand its Wick contraction in terms of the hopping parameter Counting the winding number in temporal direction, we get the fugacity expansion of the operator where z 0 is given in (3.18). The fugacity expansion coefficient is given by the expectation value Once we have two coefficients Z C (n, T, V ) and O C (n, T, V ), the VEVs of the operator with the canonical and the grand canonical ensemble are available. The canonical ensemble VEV is given by taking the ratio of two coefficient The result for the chiral condensate − d 3 x ψ ψ C /(V T 3 ) is given in the left panel of figure 10 as a function of the quark number. A VEV in the grand canonical ensemble is given by taking fugacity summation with real chemical potential  The chiral condensate − d 3 x ψ ψ G /(V T 3 ) is given in the right panel of figure 10 as a function of the real chemical potential µ/T . The condensate in the figure is a bare quantity without renormalization. Since we adopted the Wilson fermion, we have an additive correction for ψ ψ , which is not subtracted in this paper. From the figure, the chiral restoration phase transition at finite chemical potential seems to be seen. A relatively large value at around µ/T = 0 starts to decrease slowly at chemical potential µ/T ∼ 1 for low temperature. The would-be transition parameter µ c /T becomes larger for lower temperature and smaller for high temperature although we need to notice the quark mass is not the same for each β. The position of µ c /T seems to be consistent with that of the first peak for the second cumulant in figure 7 below T c indicated by arrows with corresponding color in the right panel. The finite chiral condensate above T c is an artifact due to the chiral symmetry breaking effect caused by the Wilson term and heavy quark mass. Along the vertical line at µ = 0 the lattice spacing is varied when changing the temperature so notice that the magnitude of the lattice artifact is not the same for different temperatures.
The quark number density is also given in terms of the hadronic observable d 3 ψ † ψ /(V T 3 ), which gives consistent result with that in section 5 up to renormalization factor as is plotted in figure 11.

Conclusion
We summarize what has been achieved in this paper to explore the QCD phase diagram: using the simple fugacity expansion formula derived from the hopping parameter expansion, combined with the multi-precision numerical calculation [36], we evaluate the canonical partition functions up to the large quark number. This allows us to calculate high-density region, i.e., µ/T > 1. In short, we prove that the canonical approach is the most natural way to extrapolate the imaginary chemical potential method, and indeed it is possible to study the QCD phase transition regions in a first-principle calculation.
To achieve this task, we performed the fugacity expansion of the grand partition function by using the hopping parameter expansion. This procedure seems to be valid for quark numbers around n ∼ 100 for 8 3 × 4 lattice and m π /m ρ > 0.7. One of the virtue of the hopping parameter expansion is its low computational cost. We are able to cover a very JHEP02(2016)054 wide temperature region 0.65T c ≤ T ≤ 3.62T c both in the confining and deconfining region. Note that the method is also applied to the numerator of VEVs of hadronic operators, such as meson propagators, using eq. (6.2), since hadronic operators can be expressed as the hopping parameter expansion form, whose combination with the fermion determinant has the form of the fugacity expansion. Taking summation over quark numbers we get a VEV at the real chemical potential. As an example we evaluate the pressure, the quark number density and its second cumulant and the chiral condensate. These observables show a phase transition like behavior at finite chemical potential for low temperature region, which may correspond to the confinement-deconfinement or the chiral restoration phase transition. This feature will be clarified in the near future with a high statistic simulation. We would like to add a comment that we calculate the Polyakov loop in the same procedure and observe its increase as a function of µ/T below T c , which shall be reported in the near future.
One of the biggest anxiety in this paper may be an overlap problem due to the reweighting (3.16). We performed a test of this problem by adopting µ 0 /T = 0.5 for configurations generation. We find all the results are consistent within the statistical error as we shall report in the forthcoming paper.