Study of lattice QCD at finite baryon density using the canonical approach

At finite baryon density lattice QCD first-principle calculations can not be performed due to the sign problem. In order to circumvent this problem, we use the canonical approach, which provides reliable analytical continuation from the imaginary chemical potential region to the real chemical potential region. We briefly present the canonical partition function method, describe our formulation, and show the results, obtained for two temperatures: $T/T_c = 0.93$ and $T/T_c = 0.99$ in lattice QCD with two flavors of improved Wilson fermions.


Introduction
The phase structure of hadronic matter at finite temperature and density is studied in experiments at most important modern accelerators RHIC (BNL) [1], LHC (CERN) [2] and future experiments FAIR (GSI) and NICA (JINR) will be devoted to such studies. The lattice QCD numerical simulations have a mission to provide theoretical results from the first principle calculations. Indeed, at finite temperature with zero chemical potential, the phase structure was satisfactorily investigated. But it is very difficult to study the finite density regions by the lattice QCD because of the sign problem: the fermion determinant det ∆(µ B ) at non-zero baryon chemical potential µ B is in general not real. It is impossible to apply standard Monte-Carlo techniques to computations with the partition function where S G is a gauge field action, µ q = µ B /3 is quark chemical potential, T = 1/(aN t ) defines temperature, V = (aN s ) 3 is volume, a is lattice spacing, and N t , N s -number of lattice sites in time and space directions respectively. There have been many trials, see reviews [3][4][5], and yet it is still very hard to get reliable results at µ B /T > 1. In this study we apply canonical approach. The canonical a e-mail: nikolaev.aa@dvfu.ru approach was studied in a number of papers [6][7][8][9][10][11][12][13]. We suggest new method to compute canonical partition function Z C (n, T, V), which allows to compute it for large values of n, where n is a net number of quarks and antiquarks. Our results for Z C (n, T, V) obtained with the new method are in a good agreement with results obtained with the known method of hopping parameter expansion.

Details of the new method
The canonical approach is based on the following relations. First, this is a relation between the grand canonical partition function Z GC (µ q , T, V) and the canonical one Z C (n, T, V): where ξ = e µ q /T is the fugacity and the eq. (2) is called fugacity expansion. The inverse of this equation can be presented in the following form [14]: In the right hand side of eq. (3) we see the grand canonical partition function Z GC (θ, T, V) for imaginary chemical potential µ q = iµ qI ≡ iT θ. It is known that standard Monte-Carlo simulations are possible for this partition function since the fermionic determinant is real for imaginary µ q . The QCD partition function Z GC for imaginary chemical potential is a periodic function of θ: Z GC (θ) = Z GC (θ + 2π/3). This symmetry is called Roberge-Weiss symmetry [15]. As a consequence of this periodicity the canonical partition functions Z C (n, T, V) are nonzero only for n = 3k. QCD possesses a rich phase structure at non-zero θ, which depends on the number of flavors N f and the quark mass m [16].
Quark number density n q for N f degenerate quark flavours is defined by the following equation: It can be computed numerically for imaginary chemical potential. Note, that for the imaginary chemical potential n q is also purely imaginary: n q = in qI . From eqs. (2) and (4) it follows that densities n q and n qI are related to Z C (n, T, V) (below we will use the notation Z n for the ratio Z C (n, T, V)/Z C (0, T, V)) by equations n qI /T 3 = N 2 n>0 nZ n sin(nθ) 1 + 2 n>0 Z n cos(nθ) , where N = N 3 t /N 3 s is a normalization constant. With the help of equations (5) and (6) one can construct quark number density at real and imaginary chemical potential once Z n are known.
One can compute Z GC (θ, T, V) using numerical data for n qI /T 3 via numerical integration over imaginary chemical potential: CONF12 where we omitted T and V from the grand canonical partition function notation. Then Z n can be computed as In the present work we use modified version of this approach. Instead of numerical integration in (7) we fitted n qI /T 3 to theoretically motivated functions of µ qI . It is well known that in the confining phase the hadron resonance gas model provides good description of the chemical potential dependence of thermodynamic observables [17]. Thus it is reasonable to fit the density to a Fourier expansion This type of the fit was used in Ref. [18] and conclusion was made that it works well. As an alternative one can also try to fit n qI to an odd power polynomial of θ:

Numerical results
To demonstrate our method we make simulations of the lattice QCD with N f = 2 clover improved Wilson quarks and Iwasaki improved gauge field action: with β = 6/g 2 , c 1 = −0.331, c 0 = 1 − 8c 1 , W µν denoting the plaquettes, and where P µν is the clover definition of the lattice field strength tensor and We simulate 16 3 × 4 lattices at temperatures T/T c = 0.99 and 0.93 in the confinement phase along the line of constant physics with m π /m ρ = 0.8. All parameters of the action, including c S W value, were borrowed from the WHOT-QCD collaboration paper [19]. We compute the number density on samples of N con f configurations with N con f = 1800, using every 10-th trajectory produced with Hybrid Monte-Carlo algorithm.
We employ the hopping parameter expansion to compute Z n and compare with Z n values obtained with our new method. The Wilson Dirac operator from (1) may be written in the form EPJ Web of Conferences The expansion in (16) is in fact expansion over the closed paths on the lattice, and thus it can be rewritten as where n is number of windings in the temporal direction, W n are complex coefficients which are called winding numbers. They satisfy the property W −n = W * n . In the case of the imaginary chemical potential the formula (17) will become It is important to note that hopping parameter expansion (16) converges properly only for heavy quark masses [13]. Our simulations were performed for the quark masses from this range. Below we present our results obtained below T c . In Fig. 1 we show n qI for θ ∈ [0; π/3] together with fits to eq. (9), the fit results are presented in Table 1. We found good fit with n max = 1 for T/T c = 0.93 while for T/T c = 0.99 fit with n max = 2 is necessary. In Table 1 we also show coefficients a 1 and a 3 of the Taylor expansion of eq. (10) as well as the respectve results from [19] (two last columns). We observe good agreement for the first Taylor expansion coefficient within error bars and substantially smaller error bars in our study than in [19]. For the second coefficient we see a noticeable disagreement, which means that at T < T c the fit (10) works well only up to the first order and at small enough chemical potential values.
wheref 3n = (N 3 s f 3 )/(N 3 t 3n). In the case n max = 1 this can be expressed as where I n (x) is the modified Bessel function of the first kind. Results for Z n at T/T c = 0.93 are presented in Fig. 2.
We first check if computed Z n reproduce the data for n qI /T 3 and find nice agreement between data and values of n qI /T 3 computed via eq. (6). The deviation for the full interval [0.0; π/3] is less than 0.3%. Next we compare with hopping parameter expansion, respective results are also presented in Fig. 2. We used full statistics (1800 configurations at µ qI = 0) and W n up to n = 15 in eq. (17) for this computation. One can see agreement between two results up to n = 21. We believe that the disagreement for higher n is explained by inaccuracy in computation of Z n computed by HPE. The statistical errors for Z n grow very fast with n. It is necessary to improve the HPE method accuracy before the conclusion about agreement at large n can be made.  (9) and (10). The 6th column with χ 2 and N do f demonstrates the results for the fit (9).     Table 1.
Still our result indicates that at T/T c = 0.93 the fit function (9) provides correct values of Z n and thus its analytical continuation should be valid up to values of µ q /T beyond Taylor expansion validity range. Precise determination of the range of validity of this analytical continuation will be made in future after getting more precise results for the HPE method.
Let us note that one can derive a recursion relation for Z n when n qI is presented by the function (9) with finite n max (for the derivation see Appendix). In particular for n max = 1 the recursion is just a recursion for I n which is of the formf Also one can get asymptotics for Z n at large n. For n max = 1 it is where B is some constant. This asymptotics is shown in Fig. 3 for T/T c = 0.93 with constant B obtained by the fitting over the range 400 < n < 600: B = 0.02813((1) for T/T c = 0.93. For n max = N it is different: see Appendix for the derivation. From this asymptotics it follows in particular that the coefficient f 3N has to be positive, otherwise the condition of positivity of Z 3n will not hold. Our current fitting function for T/T c = 0.99 which has f 3N ≡ f 6 < 0 does not satisfy this requirement. We need to improve statistics for this temperature to obtain the coefficient for the next harmonics in eq. (9). Evidently with present fitting function we can not go to large values of µ q for this temperature. At the end of this section we show in Fig. 4 the ratio n q /(µ q T 2 ) as a function of µ 2 q for negative and positive values. This way of presentation, borrowed from [20], allows to show in one plot the simulation results obtained at µ 2 q < 0 and analytical continuation of our fitting functions to µ 2 q > 0. For the temperature T/T c = 0.93 results seem to be reasonable, but for the temperature T/T c = 0.99 we need to improve the statistics.

Conclusions
We have presented new method to compute the canonical partition functions Z n . It is based on the fitting of the imaginary number density for all values of imaginary chemical potential to the theoretically motivated fitting function, which is Fourier-type fit (9) in the case of confinement phase.
Using fit results we compute canonical partition functions Z n ≡ Z C (n, T, V)/Z C (0, T, V) at T/T c = 0.93 and 0.99 via Fourier transformation (8). It was necessary to use the multi-precision library [21] to compute Z n which change over many orders of magnitude. For both temperatures we have checked that precision of computation of Z n was high enough to reproduce the imaginary number density n qI via eq. (6).
At the temperature T/T c = 0.93 we compared our results for Z n with Z n computed by the hopping parameter expansion. We found that two sets of Z n computed by completely independent methods agree well, see the Fig. 2. This means that the fitting function used is a proper approximations for the imaginary number density in the full range of µ qI values. Furthermore, this means that the analytical continuation to the real chemical potential can in principle be done beyond the Taylor expansion validity range since this analytical continuation of the quark number density coincides with n q computed with the help of correctly determined Z n via eq. (5). Thus the new method in principle allows to compute the number density n q beyond Taylor expansion. CONF12 Note that our new method is not limited to the heavy quark mass values like HPE, nor small µ values like Taylor expansion. Once we calculate Z n using new method, we can calculate many thermodynamical quantities, i.e. pressure, number density and its higher moments.
Using our results for the number density n qI we computed the Taylor expansion coefficients for the number density from which respective coefficients for the pressure may easily be restored. We found good agreement for the first coefficient with earlier results obtained in [19] via direct computation of these coefficients. Moreover we found that our error bars for these coefficients are in general substantially smaller than error bars quoted in [19]. Thus we confirmed analogous observation made in [20].
We obtained asymptotics of Z n at large n which is defined by the eq. (23). Such a decreasing of Z n is fast enough to provide convergence of the infinite sums in eqs. (2) and (5).