Lattice Study of QCD Phase Structure by Canonical Approach

The canonical approach is a powerful tool to circumvent sign problem in LQCD. Although it has its own difficulties it provides opportunity to determine QCD phase transition line. Using improved Wilson fermions we calculated number density at nonzero imaginary chemical potential for confinement and deconfinement phases, restored canonical partition functions Zn and did extrapolation into the real chemical potential region. We computed the higher moments of the baryon number including the kurtosis, and compared our results with information from relativistic heavy ion collision experiments.


Introduction
Many studies have tried to reveal the properties of strongly interacting quark-gluon/hadron matter from experimental and phenomenological analyses of high-energy heavy-ion collisions [1][2][3][4]. It is expected that these studies will lead to understanding of the phase diagram in the temperature -baryon density plane, which is also looked for in cosmological research. The first principle calculations, based on lattice QCD, have a potential to provide reliable fundamental information in this active area of research. However, to obtain this information, we must first overcome the "sign problem", which is described below.
The lattice QCD is a simulation study based on the grand canonical partition function, where µ is the chemical potential. Here det ∆ is the fermion determinant satisfying the relation [det ∆(µ)] * = det ∆(−µ * ).
In Monte Carlo simulations, the gluon fields, U, are generated with the probability proportional to the integrand in Eq. (1), and therefore, if det ∆ is complex, the simulations cannot be conducted. If we separate out the phase factor, i.e. rewrite the integrand as and include only the absolute value into the probability then the observables include the phase and oscillate. This makes the simulation practically impossible, and is called the "sign problem". In order to circumvent this obstacle, many approaches have been pursued, see [5] for recent review. In recent publications [6][7][8][9] where higher order cumulants were evaluated for nearly physical quark masses, mostly Taylor expansion method was employed and simulations were performed at zero chemical potential. Monte Carlo simulations for pure imaginary µ are free from the complex measure problem, as can be seen from Eq. (2). The question is how can one extract data for real µ?
The grand canonical partition function is related to the canonical partition function, Z C (n, T, V), as follows: where ξ = e µ/T is the fugacity,N is an operator of a conserved quantum number such as a baryon number or electric charge and we introduced abbreviation Z n for Z C (n, T, V). In the real simulations, we must truncate the fugacity expansion. In this letter, we are mainly concerned with the baryon number case, and we write the chemical potential µ B . For imaginary µ B (µ B = iµ I and θ I = µ I /T ), we can calculate Z n by the inverse Fourier transformation [10] as Note that Z n = n| e −Ĥ T |n ≥ 0 does not depend on µ B , and therefore one can evaluate the grand canonical partition function, Z GC , in Eq.(4) for any µ B (imaginary or real) once Z n are known.
The formula Eq.(4) is exact, and in [15], it is proved that on the finite lattice, Z GC is expressed as a finite series of the fugacity expansion Now we have a route from the imaginary to the real chemical potential regions: • Step 1: Using Eq. (5), we calculate Z n from Z GC computed at the imaginary µ B .
• Step 2: Using these Z n in Eq. (4), we construct Z GC for the real µ B .
To search for the phase transition signals one can use the moments λ m , which can be also extracted from results of the heavy ion collision experiments: Especially, λ 2 (susceptibility), λ 3 , and λ 4 , provide useful information on the phase structure. In this paper we investigate the potential of the canonical ensemble approach to reveal the QCD phases.
We can evaluate the baryon number density n B directly, for any value of the imaginary chemical potential, using standard lattice QCD algorithms: Note that the number density in imaginary chemical potential regions is pure imaginary.
On the other hand, the number density is connected with the canonical partition function, Z n as where we used Eq. (4) and relation Z n = Z −n . The direct way to extract the canonical partition functions Z n from the lattice data for n B is to fit it to Eq. (8) with Z n as fitting parameters. We tried to do it and realized that the fit goes quite unstable and some Z n 's are negative. The difficulty of fitting comes from the drastic cancellations in both the numerator and denominator in Eq. (8).
More promising way is to construct the grand canonical partition function from n B . Integrating number density over imaginary µ I , at fixed temperature T , we have where we use the fact that n B is pure imaginary and denote n BI = Im n B . We calculate Z n by inserting this Z GC into Eq (5). Then one can construct Z GC as Z GC = Z n ξ n at real µ B . This procedure provides a new method to study physics in the real chemical potential region via Monte Carlo simulations of the pure imaginary chemical potential [16,17].
There is no Ansatz until this point; therefore, Eq. (9) is exact and theoretically the calculation for any value of the chemical potential is possible. In practice, however, we must introduce some assumptions, and consequently, the reliable range of the real chemical potential values is restricted.
One way to evaluate the right hand side in Eq. (9) is to calculate the number density for many values of µ I and complete the numerical integration. In order to obtain a reliable result, we need hundreds of different µ I values, but this is computationally expensive task. In this letter we employ a simple approach -we fit the numerical data for n B and use the fit function in Eq. (9). Thus the idea is to fit the number density as a function of µ I to Ansatz Eq. (10) and (11) presented below, compute the partition function for imaginary µ B from Eq.(9), then compute Z n using Eq. (5) This idea is a continuation of another concept Refs. [18][19][20] -fit the data at imaginary µ B and do analytical continuation to the real axis.
The authors of Ref. [8] have recently reported a thorough analysis. In Refs. [18][19][20], the authors pointed out that the number density for the imaginary chemical potential is well approximated by a Fourier series at T < T c , and by a polynomial series at T > T c 1 , In Refs. [16,17] we confirmed these conclusions with higher precision. Therefore in the present study in computation of Z GC (θ)/Z GC (0) we use parameterizations Eqs. (10) and (11). We obtained the best fits for k max = 2 with f 3 = 0.0871(3) and f 6 = −0.00032(27) (χ 2 /do f = 0.93) at T/T c = 0.93, a 1 = 1.5570 (7), a 3 = −0.3300(13) (χ 2 /do f = 0.67) at T/T c =1.35 for our data n BI . We also studied systematic error due to parametrization function -for confinement phase we parametrized the number density not only with Fourier anzats but also with polynomial and varied k max .  To estimate the statistical error, we apply a version of a bootstrap method: Here one bootstrap sample consists of a set of standard bootstrap samples of the number density created for every value of µ I . On each bootstrap sample we estimate n B with Eq. (10) (T < T c ) or Eq. (11) (T > T c ) (coefficients f 3k or a 2k−1 are different on each sample), and calculate Z n according to Eqs. (9) and (5). Then, using these Z n and Eq. (6) with number of terms restricted to n max , we calculate observables.
It is important to note that the baryon number density fitted to Eqs. (10) or (11) can be analytically continued to the real chemical potential values -it was previously studied in Refs. [16,20]. Here we also repeat this procedure to illustrate drawback of this approach.
As it is shown at figures 1 and 2 the number density can be fitted to different functions at imaginary µ but when we extrapolate them to real region they give substantially different results for high enough µ. There is no strict argument how to choose correct parametrization function so reliable µ range (where all fitting function predict same results) can be small. As one can see from Fig.2 the systematic error due to uncertainty in the choice of the fitting function is small for mu B /T < 3. But there might be other suitable fitting functions which introduce even larger systematic error for real mu.
In opposite, the canonical approach provides us with useful information on the range of reliability at real mu. Because of Z n =< n|e −Ĥ T |n >, Z n must be positive. If Z n become negative for some n > n max it means that our Ansatz used for Z n calculation can not describe physics for these n. In the statistical analysis as described above we used only positive Zn. It is important to note that in each Bootstrap sample n max can be different due to slightly different coefficients in Ansatz.
In the Fig. 3 and 4 one can see Z n /Z 0 calculated with different number of terms k max in Eq. (10) and (11). First of all, for small n all parametrization functions give same Z n but for higher n they deviate. For polynomial fits starting from some large n = n max value the sign of Z n starts to alternate while the absolute value changes very slowly. We think these fits should be rejected since Z n should be positive. For fits (10) with even k max similar thing happens: starting from some n max there appear alternation of the sign of Z n . The difference from the polynomial fit is that the sign alternate not at   (11)) with different number of termsk max . Temperature T/T c = 0.93 (Confinement). In the legend number is equal to k max , 'sin' denotes Fourier type fits Eq. (10), 'pol' denotes polynomial fits Eq. (11)) every step in n as seen in Fig. 3 for k max =10. We thus conclude that the Fourier fits with even k max should be also avoided or used for the range of n values below respective n max .
It is important to note that all polynomial functions give small n max ≈ 50 despite the fact that they also give good fit of our data. We see from Figs. 3 and 4 that dfferent number of terms in Fourier series give rise to different Z n which agree within error bars. The inceasing of number of k max increases errors of coefficients in (10) and as a result errors of Z n are also increased. We think the natural choice for this temperature is k max = 2. In addition Z n analisys agreeds with known fact that in the confinement phase baryon density is well described by Hadron Resonance Gas model and therefore is approximated by sum of sines.

Moments of the baryon number density
In Ref. [21], the canonical partition functions, Z n , were extracted from the RHIC experiments data. With these Z n , we construct moments of the baryon number density for different RHIC energy values and compare with our results.
In the relativistic heavy-ion collision experiments, λ 2 /λ 1 is expected to be good indicators for detecting the QCD phase transition [2,22]. In Fig. 5 we show the ratio λ 2 /λ 1 as calculated by the integration method described above together with those extracted from the RHIC STAR data.  We see from the Fig.5 that at T > T c our results deviate from RHIC data, while the confinement phase results can be consistent with the experimental data with exception for √ s NN =11.5. Estimated temperature at √ s NN =11.5 in Ref. [23] is significantly below other energy data.
Note that in Ref. [21] Z n were constructed from the proton multiplicity data, not the baryon multiplicity. Therefore, the results should be considered as a proxy for the real baryon number moments. Nevertheless, in the confinement regions we see very good agreement for λ 2 /λ 1 between the lattice calculation and those estimated from RHIC data. It is important since the fact that for RHIC data we have only small Z n number, i.e., n max reliability range for it is much smaller, we estimatre it as µ reliable /T ∼ 1.2 − 2.5 depending on √ s NN .

Concluding Remarks
In this letter, we study an approach for revealing the QCD phase structure using lattice QCD simulations. Prior to this study, it was believed that this was impossible because of the sign problem; only small density regions could be studied by extrapolating from the data at µ B = 0. However, all relevant information on the QCD phase at finite baryon density is contained in the imaginary chemical potential regions, 0 ≤ µ I /T ≤ π. The question is how to map this information to the real chemical potential. Eq. (4) provides a possible solution, because Z n can be calculated in the imaginary chemical potential regions. Since numerical Monte Carlo simulations provide results with finite accuracy, we should find practical methods which work. We parametrize the number density in the imaginary chemical potential with the Ansatz (the Fourier series in the confinement and polynomial series in the deconfinement) and integrate them to get the grand partition function. Z n and other observables are then calculated from them.
This method produces Z n up to n max which is determined by used Ansatz and current statistics. However, this is not the first principle calculation due to introducing an assumption to the number density, analyzing n max behavior we can estimate reliability of this assumption.
On the other hand, better interpolation procedure (cubic spline for example) or numerical integration rather than any Ansatz release our data from assumptions. This study will be reported in future.
We studied how errors of calculations made at the imaginary chemical potential propagate into errors at the real one. We also found range of reliability of our results by imposing the condition that Z n have to be positive.
We then investigate whether we can estimate λ 2 /λ 1 . The results are consistent with the values estimated from the RHIC experiments as shown in Figs. 5 and ??. This is very encouraging. More realistic simulations with the physical quark mass and small lattice spacing, may even predict the temperature of the experimental data.
When we map the information from pure imaginary to real chemical potential by the canonical method, the reliable regions for baryon density is µ reliable /T ∼ 4 − 5 at T/T c = 0.93 and µ reliable /T ∼ 2 − 2.5 at T/T c = 1.35. This limitation comes from the finite numbers of Z n as a consequence of polynomial (in the deconfinement phase) /Fourier (in the confinement phase) anzats. They can be increased by increasing statistics.