Lattice QCD at finite baryon density using analytic continuation

We simulate lattice QCD with two flavors of Wilson fermions at imaginary baryon chemical potential. Results for the baryon number density computed in the confining and deconfining phases at imaginary baryon chemical potential are used to determine the baryon number density and higher cumulants at the real chemical potential via analytical continuation.


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 baryon density -temperature plane. New experiments will be carried out at FAIR (GSI) and NICA (JINR). To explore the phase diagram theoretically it is necessary to make computations in QCD at finite temperature and finite baryon chemical potential. For finite temperature and zero chemical potential lattice QCD is the only ab-initio method available and many results had been obtained. However, for finite baryon density lattice QCD faces the so-called complex action problem (or sign problem). Various proposals exist 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.
Here we consider the analytical continuation from imaginary chemical potential.
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. It is known that the standard Monte Carlo simulations are possible for the grand canonical partition function Z GC (θ, T, V) for imaginary chemical potential µ q = iµ qI ≡ iT θ. since the 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 [6]. 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 temperature 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. [7]. 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 . In this work 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 at temperature T > T RW . This type of the fit was also used in Refs. [8][9][10][11].
In the confining phase (below T c ) the hadron resonance gas model provides good description of the chemical potential dependence of thermodynamic observables [12]. Thus it is reasonable to fit the density to a Fourier expansion Again this type of the fit was used in Refs. [9,10] and conclusion was made that it works well. We use both types of the fitting function in the deconfining phase at T c < T < T RW . We made simulations of the lattice QCD with N f = 2 clover improved Wilson quarks and Iwasaki improved gauge field action. The more detailed definition of the lattice action can be found in Ref. [8]. The simulations were made on 16 3 × 4 lattices. We obtained results at temperatures T/T c = 1.35, 1.20, 1.08, and 1.035 in the deconfinement phase and 0.99, 0.93, 0.84 in the confinement phase along the line of constant physics with m π /m ρ = 0.8. We also present here our preliminary results for smaller quark mass with m π /m ρ = 0.65 At this quark mass the simulations were made at T/T c = 1.32, 1.18, 1.07, 1.00, 0.94, 0.86. The parameters of the action, including c S W value were borrowed from the WHOT-QCD collaboration paper [13]. We compute the number density on samples of N con f configurations with N con f = 1800 or 3800, using every 10-th trajectory produced with Hybrid Monte Carlo algorithm.

Quark number density
In this section we compare results for the quark number density obtained at the imaginary chemical potential for two values of the quark mass. In Figure 2 the data for the deconfinement phase are shown. One can see that at small values of µ qI /T the number density for two quark masses differ only slightly for comparable values T/T c . At the same time effects of the quark mass decreasing are quite visible in the range µ qI /T > 0.8. As can be seen from Figure 3 in the confinement phase the differences between results for two quark masses might be mostly due to differences in the T/T c values.
This assumption is supported by comparison of the virial coefficients f n depicted in Figure 4. Note logarithmic scale for Y-axes in this Figure.

Taylor expansion coefficients
The Taylor expansion coefficients for the pressure are introduced as follows: where ∆P(T, µ B ) = P(T, µ B ) − P(T, 0), P 2k are Taylor expansion coefficients, χ 2k are called generalized susceptibilities. Before we discuss the Taylor expansion coefficients few comments about the fits follow. The fits to function eq. (5) are very stable with respect to change of the fitting range. This was observed for T/T c = 0.84, 0.93. Since the statistical error for the coefficient f 3 is small at low temperatures we obtain the Taylor coefficients a k with low error even for high values of k. We should note that there is a source of uncertainty which we cannot estimate reliably. This is the contribution from the higher terms in the Fourier decomposition eq. (5). Our data indicate that for low temperatures f 6 is at least factor 100 smaller than f 3 . This implies that for Taylor coefficients a 1 and a 3 the contribution from the second term in eq. (5) should be small while starting from a 5 this contribution might be substantial. For consistency check of our results we also applied polynomial fit at low temperatures. The polynomial fit was applied for restricted range of µ qI values. We obtained results compatible with respective Taylor coefficients a k , k = 1, 3, 5 within error bars for T/T c = 0.84, 0.93. For T/T c = 1.035 and 1.08 the results for a k obtained with two kind of fits are in agreement for k = 1 only. We then used a 3 and a 5 values obtained from the polynomial fit. For the lighter quark mass we made similar computations.
The generalized susceptibilities χ n for n = 2, 4, 6 which are proportional to the Taylor coefficients P n are presented in Figure 5 for m π /m ρ = 0.8 and in Figure 6 for m π /m ρ = 0.65. We should note that our results for P 2 and P 4 are in agreement within error bars with results of Ref. [13] where direct computation of the Taylor coefficients was done. But our error bars are substantially lower.
The Taylor coefficients were recently computed for the physical quark masses on lattices with small lattice spacing (and even in the continuum limit) in Refs. [11] and [14]. In Ref. [11] the analytical continuation was used while in Ref [14] direct method was employed. Results of these two computations were found to be in a good agreement [14]. Our results presented in Figure 5 and Figure 6 are in good qualitative agreement with results of Refs. [11,14] for all three susceptibilities. Quantitatively our results at T > T c are substantially higher what should be expected from the large lattice spacing effects estimated in Ref. [13].  In Figure 7 and Figure 8 we show the ratios χ B 4 /χ B 2 and χ B 6 /χ B 2 , respectively. One can see that for these ratios our results are in very good agreement with results of Ref [14] taken from their Figure 3. Substantial difference is observed only for χ B 4 /χ B 2 at 1 < T/T c < 1.1. This agreement indicates that the finite lattice spacing effects are substantially cancelled in the ratios of the susceptibilities.

Conclusions
We computed the baryon number density and generalized susceptibilities χ n in the lattice QCD with two flavors of Wilson fermions on 16 3 × 4 lattices using analytical continuation. Comparing results for two quark masses we found that they do not differ substantially. Comparison of our results for the ratios χ B 4 /χ B 2 and χ B 6 /χ B 2 with results of Ref [14] where simulations were done at the physical quark masses and small lattice spacing we found surprising agreement. This agreement indicates that our recent results [15] showing agreement of cumulant ratios computed on the lattice with respective experimental results are not accidental.