Study of lattice QCD at finite chemical potential using canonical ensemble approach

New approach to computation of canonical partition functions in $N_f=2$ lattice QCD is presented. We compare results obtained by new method with results obtained by known method of hopping parameter expansion. We observe agreement between two methods indicating validity of the new method. We use results for the number density obtained in the confining and deconfining phases at imaginary chemical potential to determine the phase transition line at real chemical potential.


Introduction
Recent results of heavy ion collision experiments at RHIC [1] and LHC [2] shed some light on properties of the quark gluon plasma and the position of the transition line in the temperature -baryon density plane. New experiments will be carried out at FAIR (GSI) and NICA (JINR). To fully explore the phase diagram theoretically it is necessary to make computations at finite temperature and finite baryon chemical potential. For finite temperature lattice QCD is the only ab-initio method available and many results had been obtained. However, for finite baryon density lattice QCD faces the socalled complex action problem (or sign problem). Various proposals exist at the moment to solve this problem see, e.g. reviews [3][4][5] and yet it is still very hard to get reliable results at µ B /T > 1. Our work is devoted to developing the canonical partition function approach.
The fermion determinant at nonzero baryon chemical potential µ B , det ∆(µ B ), is in general not real. This makes 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 ) is temperature, V = (aN s ) 3 is volume, a is lattice spacing, N t , N s -number of lattice sites in time and space directions. a e-mail: vitaly.bornyakov@ihep.ru The canonical approach was studied in a number of papers [6][7][8][9][10][11][12][13]. It is based on the following relations. Relation between grand canonical partition function Z GC (µ q , T, V) and the canonical one Z C (n, T, V) called fugacity expansion: where ξ = e µ q /T is the fugacity. The inverse of this equation can be presented in the following form [14] Z C (n, Z GC (θ, T, V) is the grand canonical partition function for imaginary chemical potential µ q = iµ qI ≡ iT θ. Standard Monte Carlo simulations are possible for this partition function since fermionic determinant is real for imaginary µ q . The QCD partition function Z GC 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 nonzero θ, which depends on the number of flavors N f and the quark mass m. This phase structure is shown in Fig. 1. T c is the confinement/deconfinement crossover point at zero chemical potential. The line (T ≥ T RW , µ I /T = π/3) indicates the first order phase transition. On the curve between T c and T RW , the transition is expected to change from the crossover to the first order for small and large quark masses, see e.g. [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 q /T 3 = N 2 n>0 nZ n sinh(nθ) 1 + 2 n>0 Z n cosh(nθ) , n qI /T 3 = N 2 n>0 nZ n sin(nθ) 1 + 2 n>0 Z n cos(nθ) , where N is a normalization constant, N = Our suggestion is to compute Z n using equation (5) for n qI .
One can compute Z GC (θ, T, V) using numerical data for n qI /T 3 via numerical integration where we omitted T and V from the grand canonical partition function notation. Then Z n can be computed as In our work we use modified version of this approach [17]. Instead of numerical integration in (6) we fitted n qI /T 3 to theoretically motivated functions of µ qI . It is known that the density of noninteracting quark gas is described by We thus fit the data for n qI to an odd power polynomial of θ in the deconfining phase. This type of the fit was also used in Ref. [18] and Ref. [19].
In the confining phase (below T c ) the hadron resonance gas model provides good description of the chemical potential dependence of thermodynamic observables [20]. Thus it is reasonable to fit the density to a Fourier expansion Again this type of the fit was used in Ref. [18] and conclusion was made that it works well.
To demonstrate our method we made simulations of the lattice QCD with N f = 2 clover improved Wilson quarks and Iwasaki improved gauge field action, for detailed definition of the lattice action see [17]. We simulate 16 3 × 4 lattices at temperatures T/T c = 1.35, 1.20 and 1.08 in the deconfinemnt phase and 0.99, 0.93, 0.84 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 [21]. 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.

Results
In Fig. 2 and Fig. 3 we show numerical results for n qI as a function of θ in deconfining and confining phases, respectively [17].
Note, that behavior of n qI at T = 1.08 is different from that at higher temperatures. This temperature is below T RW and at θ = π/3 there is no first order phase transition, n qI is continuous. Instead there is a crossover to the confinement phase at about θ = 0.92(2). It is not yet clear how to fit the data over the range of µ qI covering both deconfining and confining phase. Here we use fit to function (9) with n max = 3 over the range [0, 0.8], i.e. including only the deconfining phase. In this case we should consider the fit as a Taylor expansion.
We computed Z n using new procedure described in the previous section and compared with results obtained with use of the hopping parameter expansion. We found good agreement between two methods indicating that the new method works well [17]. This allows us to make analytical continuation to the real values of µ q beyond the Taylor expansion validity range (for all temperatures apart from T/T c = 1.08). Respective results are shown in Fig. 4. Using the results for the number   density n q we can compute the temperature of the transition from the hadron phase to quark-gluon plasma phase using the following procedure. Our results for the number density at temperatures T > T c as function of the chemical potential µ q are reliable even for large µ q for T > T RW , while for T c < T < T RW they are reliable at the moment for small µ q only. We can use these results to compute the pressure ∆P decon f (T, µ q ) = P(T, µ q ) − P(T, 0) as function of µ q and then extrapolate pressure for fixed µ q to temperatures T < T c . To get good extrapolation we need pressure computed for more values of temperature than is available now but three values which we have in this work is enough to demonstrate the idea. We then find the transition temperature T c (µ q ) solving numerically equation ∆P decon f (T, µ q ) = ∆P con f (T, µ q ), where ∆P con f (T, µ q ) is pressure computed from results for the number density n q we obtained in the confinement phase. In this paper we use extrapolation for the coefficients in (9) rather than for pressure itself. This extrapolation is shown in Fig. 5.
The result for the fit parameter is C = 0.07. We have to emphasize that in Figs. 6 and 7 we do not show statistical or systematic errors since we are not yet able to compute them. These figures are shown to illustrate our idea about the method to compute the transition temperature dependence on the chemical potential. Still it is encouraging that the value of the coefficient C obtained is of correct order of magnitude. We can compare with Refs. [22,23] where C = 0.051(3) and C = 0.065 (7) respectively were found in lattice QCD with N f = 2. In both papers analytical continuation of T c (µ q ) T c from imaginary µ q was used. Thus we presented new method to compute the canonical partition functions Z n . It is based on fitting of the imaginary number density for all values of imaginary chemical potential to the theoretically motivated fitting functions: polynomial fit (9) in the deconfinement phase for T above T RW and Fourier-type fit (10) in the confinement phase. We also explained how the transition line at real µ q can be computed using results we obtained for the number density n q . We are planning to increase the number of temperature values to improve the precision of this method. It is worth to note that the method of direct analytical continuation for T c (µ q ) T c from imaginary chemical potential can be applied to our data. This will help us to improve the method presented here.