Study of QCD critical point using the canonical ensemble method
Abstract
The existence of the QCD critical point at nonzero baryon density is not only of great interest for experimental physics but also a challenge for the theory. Any hint of the existence of the first order phase transition and, particularly, its critical point will be valuable towards a full understanding of the QCD phase diagram. We use lattice simulation based on the canonical ensemble method to explore the finite baryon density and finite temperature region and look for the QCD critical point. As a benchmark, we run simulations for the four degenerate flavor QCD where we observe a clear signal of the expected first order phase transition. In the two flavor case, we do not see any signal for temperatures as low as . Although our real world contains two light quarks and one heavier quark, three degenerate flavor case shares a lot of similar phase structures as the QCD. We scan the phase diagram using clover fermions with on lattices. The baryon chemical potential is measured as we increase the baryon number and we see the characteristic “Sshape” that signals the first order phase transition. We determine the phase boundaries by Maxwell construction and report our preliminary results for the location of critical point for the present lattice.
Study of QCD critical point using the canonical ensemble method
Anyi Li^{†}^{†}thanks: Speaker. ^{†}^{†}thanks: Present address: Department of Physics, Duke University, Durham, NC 27708
Department of Physics and Astronomy, University of Kentucky, Lexington, KY 40506
Email: anyili@phy.duke.edu
\abstract@cs
1 Introduction
Quantum Chromodynamics (QCD) is a fundamental theory which describes the strong interaction of quarks and gluons—the basic blocks of matter. The knowledge of QCD at finite temperature and finite baryon chemical potential is essential to the understanding of a variety of phenomena. Exploring the QCD phase diagram, including identification of different phases and the determination of the phase transition line is currently one of the most studied topics [1, 2]. Guided by phenomenology and experiments, a candidate QCD phase diagram was proposed [3] – Fig. 1 shows a sketch of the conjectured phase diagram. In certain limits, in particular for large temperatures or large baryonchemical potential , thermodynamics is dominated by shortdistance QCD dynamics and the theory can be studied analytically. However, the most interesting phenomena lie in regions where nonperturbative features of the theory dominate. The only known systematic approach is the firstprinciple calculation via lattice QCD using Monte Carlo simulations.
The thermodynamics of strongly interacting matter has been studied extensively in lattice calculations at vanishing baryon chemical potential. Current lattice calculations strongly suggest that the transition from the low temperature hadronic phase to the high temperature phase is a continuous but rapid crossover happening in a narrow temperature interval around [4, 5, 6, 7, 8, 9, 10]. On the other side of the phase diagram — large baryon chemical potential but very low temperature, a number of different model approaches [11, 12, 13, 14, 15, 16, 17, 18] suggest that the transition in this region is strongly first order, although this argument is less robust. This first order phase transition is expected to become less pronounced as we lower the chemical potential and it should terminate at a second order phase transition point – the critical point.
The search for the QCD critical point has attracted considerable theoretical and experimental attention recently. The possible existence of this point was introduced sometime ago [13, 14]. It is apparent that the location of the critical point is a key to the understanding of the QCD phase diagram. For experimental search of the critical point, it has been proposed to use heavy ion collisions at (RHIC) [19, 20]. The appearance of this point is closely related to hadronic fluctuations which may be examined by an eventbyevent analysis of experimental results. The upcoming RHIC lower energy scan and future LHC and FAIR (GSI) will focus on the region and where the critical point is predicted to exist in theoretical models. However, in order to extract unambiguous signals for the QCD critical point, quantitative calculations from firstprinciple lattice QCD are indispensable.
Remarkable progress in lattice simulations at zero baryon chemical potential has been made in recent years; however, simulations at nonzero chemical potential are difficult due to the complex nature of the fermionic determinant at nonzero chemical potential. The phase fluctuations produce the notorious “sign problem”. The majority of current simulations are focusing on the small chemical potential region where the “sign problem” appears to be controllable. Most of them are based on the grand canonical ensemble (, as parameters). Since the existence and location of the critical point is still unknown, we need an algorithm which can be applied beyond the small chemical potential region. This is one of the motivations of our studies via the canonical ensemble approach [21, 22, 23, 24, 25].
We have proposed an algorithm based on the canonical partition function to alleviate the overlap and fluctuation problems [26, 27]. The method we use is computationally expensive since every update involves the evaluation of the fermionic determinant; however, finite baryon density simulation based on this method has been shown to be feasible [27, 28]. In this approach, we measure the baryon chemical potential to detect the phase transition. With the aid of the winding number expansion technique [29, 30, 31, 32], a program was outlined to scan the QCD phase diagram in an effort to look for the critical point [33, 34].
In this study, we present results based on simulations using lattices with clover fermions. We have run simulations using two, three and four degenerate flavors of quarks. As a benchmark, we first study the four flavor case. This system has a first order phase transition that extends all the way to zero density and the low density region has been reliable studied using staggered fermions. We map out the phase transition line using our method and compare it with the results of previous studies; we find very good agreement. For the two flavor case we present results for three different temperatures, the lowest one being . In this case we do not see any clear signal for a phase transition. While surprising, this is consistent with results from other studies [35]. The most interesting results we report are for the case, where the evidence of a critical point would not only suggest the existence of such a point in the real world but also give an estimation of its location. We scan the phase diagram and we find signals for a first order phase transition for temperatures below . We locate the critical point at the intersection of the boundaries of the coexistence region.
2 Canonical partition function
One way to show how to build the canonical ensemble in lattice QCD is to start from the fugacity expansion,
(2.0) 
where is the net number of quarks (number of quarks minus number of antiquarks) and is the canonical partition function. Using the fugacity expansion, it is easy to see that we can write the canonical partition function as a Fourier transform of the grand canonical partition function,
(2.0) 
after integrating over the imaginary chemical potential.
For our study, we use the Iwasaki improved gauge action and clover fermions [36, 37, 38].The imaginary chemical potential is introduced via a phase added to the time links from the last time slice.
As an illustration, we will consider the case of two degenerate flavors. After integrating out the fermion part, we get
(2.0) 
where
(2.0) 
is the projected determinant with the fixed net quark number . Using charge transformation and hermiticity property of the action, we can prove that
(2.0) 
This property allows us to rewrite the partition function as
(2.0) 
The integrand is real but not necessarily positive. To compute observables based on this partition function we rewrite it as:
(2.0)  
where
(2.0) 
and
(2.0) 
We separate the phase and generate an ensemble distributed according to the positive weight ; the phase factor is reintroduced in the observable. To generate this ensemble a candidate configuration is “proposed” by the standard Hybrid Monte Carlo algorithm [39, 40, 41, 42], then an accept/reject step is used to correct the weight. We note that the accept/reject step is based on the determinant ratio which alleviates the fluctuation problem [43, 23, 26] and enhances the acceptance rate [27].
3 Winding number expansion
Most of the simulation time is spent on the accept/reject step, specifically on computing determinant of the fermion matrix. On the lattice, the matrix has entries. Although it’s very sparse, exact determinant calculation is numerically demanding even on this small lattice. Moreover, the determinant needs to be evaluated many times at every accept/reject step since we need to compute its Fourier transform; our original approach was to use an approximation where the continuous Fourier transformation is replaced by a discrete one, i.e.:
(3.0) 
It was shown that the errors introduced by this approximation are small [34] for small quark numbers. However, there are two problems with this approach: (i) the computation time increases linearly with the net quark number, (ii) although evaluating the determinant times should theoretically allow us to compute projection numbers as large as , numerically we found that for large quark numbers, the Fourier components become too small to be evaluated with enough precision, even using double precision floating point numbers. It has been shown [29] that the the results of the projected determinant for larger than would differ significantly for different choice of , which signals a numerical instability. To see this problem, we use to evaluate the fermion determinant and calculate Fourier projection using discrete Fourier transform. One would expect that should decrease exponentially. We see that this is indeed the case for discrete Fourier transform at small quark number, as shown in Fig. 2 (left panel). However, as the quark number gets close to 30, the projected determinant calculated using the discrete Fourier transform approximation flattens out and stops falling below the double precision limit of . This is the onset of the numerical instability.
This happens because oscillates rapidly at larger , leading to large cancellation in the sum over . The accumulation of roundoff errors makes it impossible to evaluate the projected determinant. This numerical challenge led us to develop a new method which should be free of the above numerical problem for quark numbers relevant to the phase transition region in this study [29]. The idea of our new method is to first consider the Fourier transform of instead of the direct Fourier transform of . Using an approximation based on the first few components of , we can then analytically compute the projected determinant. The efficacy of the method can be traced to the fact that the Fourier components of are the number of terms in the expansion which characterizes the number of quark loops wrapping around the time boundary. They are exponentially smaller with increasing winding numbers. This is why we can approximate the exponent of the determinant very accurately with a few terms which, in turn, allows us to evaluate the Fourier components of the determinant precisely.
To see how this works, we look at the hopping expansion of the . We start by writing the determinant in terms of the trace log of the quark matrix
(3.0) 
It is well known that corresponds to a sum of quark loops. We separate all these loops in classes depending on how many times they wrap around the lattice in the temporal direction. We have then
where is the winding number of quark loops wrapping around the time direction and is the weight associated with all these loops with winding number . is the contribution with zero winding number. Eq. (3) can be rewritten as
where and are independent of . Using Eq. (3.0) and Eq. (3) we get
(3.0) 
The Fourier transform of the first order in the expansion can now be computed analytically; we have
(3.0) 
where is the Bessel function of the first kind with rank .
For higher orders in the winding number expansion, we compute the Fourier transform based on the truncated Taylor expansion:
(3.0)  
The projected determinant is written in terms of the linear combination of Bessel functions, the coefficients can be easily computed analytically. Using Eq. (3.0) and the recursion relation for the Bessel function, , the winding number expansion method (WNEM) can be extended to higher orders.
As we mentioned before, the efficacy of the method rests on the fact that we can get a very good approximation of the exponent using only a few terms in the Fourier expansion. However, the evaluation in Eq. (3.0) is based on a mixture of analytical integration and a truncated Taylor expansion. Since the error of the truncated Taylor expansion is not well controlled, to produce our final results, we decided to employ a method which can evaluate Eq. (3.0) more precisely and thus correct the necessary errors introduced by the Taylor approximation. The methodology can be viewed as a form of reweighting and the approximation level is then tuned to produce errors that are smaller than the statistical errors of the final result. The strategy is the following:

During the procedure of generating ensembles, we compute the determinant for 16 phases and we keep only 6 winding loops in the expansion which has been shown to be precisely enough for the approximation of [29]. The evaluation of the projected determinant is done by truncating the second and third winding loops terms to their tenth order, and first order for the remaining terms. We denote the approximated determinant using truncated Taylor series in Eq. (3.0) by .

Once the ensembles are generated, we evaluate the determinant of each configuration for phases which admits projection up to support as large us winding loops reliably. Instead of evaluating the integral Eq (3.0) by the truncated Taylor series, we evaluate it numerically through a multiprecision library (GMP library) by setting precision as 512 digits after decimal point. In order to choose , we include the higher order of winding loops until the projected determinant converges, the minimal number of winding loops included will be equal to .

Since the ensemble generated in the Monte Carlo process involves an approximation of the projected determinant and is not the target ensemble, we correct for this in the observable with the following reweighting
(3.0) where denotes the average over the ensemble generated with the approximated determinant and
(3.0) is the reweighting phase where is from the exact Fourier transform without Taylor expansion. If the approximation due to Taylor expansion is very good then will be very close to 1; if the difference between the truncated determinant and the exact one is much smaller than the error bars on then the reweighting step is not really needed and we can consider the approximation exact. On the other hand, if the average is significantly different from 1 we need to do reweighting. The worst case is when the approximation is so bad as to introduce a sign problem ( overlaps with zero within the error bars). In that case the approximation is invalid since the generated ensemble cannot be used to compute any observable.
The detailed discussion of this methodology will be presented somewhere [44]. By employing the above strategy, even if the approximation is not perfect, it can be corrected by a reweighting after the ensemble is generated. The only possible issue is to make sure that the reweighting phase doesn’t introduce too much noise in the observable; if that is the case, the approximation needs to be improved and a new ensemble needs to be generated.
Before we conclude this section, we would like to compare the merits of WNEM with the those of discrete Fourier transform. We compute the values of the projected determinant determined using the discrete Fourier transform with only in Eq. (3.0). In Fig. 2 (right panel), these results are labeled discrete (16). We see that this approximation is only valid up to . This is to be compared to WNEM (16) which takes the same computational time and yet does not suffer from this problem and, as expected, the evaluated determinant continues to decrease as the quark number is increased. It is this projected determinant evaluated from WNEM that allows us to scan a wide range of densities on the QCD phase diagram.
4 Results
4.1 Baryon chemical potential
Before starting the discussion on the QCD phase diagram, we will first present the “tool” we used to determine the phase transition in the canonical ensemble: the baryon chemical potential. We measure the chemical potential and plot it as a function of the net quark number . Due to the contribution of the surface tension to the free energy at finite volume, the chemical potential in the mixed phase region has an “Sshape” structure. We scan the phase diagram looking for this Sshape a structure to signal a first order phase transition.
In the canonical ensemble, the baryon chemical potential is measured by taking the difference of the free energy after adding one baryon, i.e.
(4.0) 
where
(4.0)  
(4.0) 
are phase factors measured in the ensemble with baryon number where in Eq. (4.0) stands for the average over the ensemble generated with the measure .
4.2 Phase diagram for and
We turn now our attention to the QCD phase diagram. Although it has been speculated for quite a some time that the QCD phase diagram may look like Fig 1, unfortunately, there is little solid evidence to to support this scenario. The finite temperature studies at zero baryon density have shown that the transition from hadronic matter to quark gluon plasma is a rapid crossover. Lattice QCD evidence regarding the existence and location of the second order transition point will help tremendously in sharpening our understanding of the QCD phase diagram.
Before we present our results, we would like to point out the difference between the phase diagram in the grand canonical ensemble and the one in the canonical ensemble. We plot the canonical partition function phase diagram in Fig. 3. In terms of and , the first order phase transition line becomes a phase coexistence region which has two boundaries that separates it from the pure phases. As we approach the critical point, the two boundaries will eventually meet at that point. We determine the boundaries of the coexistence phase at different temperatures and locate the critical point by extrapolation.
As a first attempt, we carry out simulations with four and two degenerated quark flavors. The reason we start with two and four flavors is the following: first of all, it is easier to simulate an even number of flavors with HMC [40]. Secondly, in the fourflavor case, the first order phase transition line extends all the way to zero density. This provides a testing ground for the algorithm since there are reliable results mapping the transition line in the small chemical potential region. Thirdly, the twoflavor case is expected to have a phase diagram similar to that of real QCD. The expected phase diagrams for two and four flavor cases are shown in Fig. 3.
For this reason, we will use the four flavor simulation as a benchmark to demonstrate the methodology of determining the phase boundaries and locating the critical point before tackling the more realistic case. The expected phase diagrams for two and four flavor cases are shown in Fig. 3.
Given the first order transition at zero chemical potential for , it is natural to ask whether this first order is preserved when we switch on the chemical potential. As we mentioned before, the first order phase transition is reflected in an “Sshape” structure in the vs. plot (or vs. plot as we will present). To proceed in this way, we decide to scan the phase diagram by fixing the temperature below while varying .
We scan the phase diagram by fixing the temperature below while varying . Once we cross into the coexistence region, in finite volume, the nonzero contribution from the surface tension causes the appearance of a “doublewell” effective potential [35] which leads to an Sshaped behavior in the chemical potential plot. In the thermodynamic limit, the surface tension contribution goes away since it is a surface term and the free energy scales with the volume; the chemical potential will then stay constant in the mixedphase region.
We present our results for in Fig. 4 for three different temperature below . To identify the boundaries of the coexistence region and the critical value for the baryon chemical potential, we rely on the “Maxwell construction” [25]: we select several points in the Sshape region and fit them with a thirdorder polynomial. We determine the critical chemical potential and the boundary points and by requiring the two areas between the fitted curve and a constant horizontal line which intersects the curve at and to be equal. For all reasonable fits, we find the values of the boundary points, and , and the value of the critical chemical potential () to be fairly insensitive to our choice of the fit function or fit region. The simple third order polynomial fit is sufficient.
We perform the Maxwell constructions for the three temperatures we studied: , and ; the results are presented in Fig. 4. Having determined the and for three different temperatures, we plot boundaries of the coexistence region by a simple extrapolation. Although there is no "critical point" for the case, it is expected that the two coexistence phase boundaries should cross at zero chemical potential and . To determine the crossing point, we fit the boundary lines using a simple even polynomial in baryon density
(4.0) 
to do the extrapolation. The phase boundaries and their extrapolation are plotted in Fig. 5. We find the intersection point at and which is consistent with our expectation: and .
We would like to point out that the critical temperature was determined using a different set of simulations. We run simulations at zero baryon density and we monitored the Polyakov loop susceptibility as we increased the temperature. These simulations were performed using the standard HMC algorithm and we used two lattice sizes: and . The critical temperature was determined using the peak of the Polyakov loop susceptibility and we also checked that the susceptibility peak scales with the volume as expected. The fact that the intersection point is close to the critical temperature determined from a set of simulations using a different methodology constitutes a convincing crosscheck of the correctness of our method.
We compare our results to the study of analytic continuation from the imaginary chemical potential with stagger fermions [45]. The results are presented in Fig. 6 (left panel) – we see that the agreement is quite good in spite of the fact that the quark mass used in the staggered fermion study, , is significantly different from the one we use. This is not surprising since the phase transition curve for is fairly independent of the quark mass [46]. As a check that the reweighting does not change drastically the observable results, we also plot the results with no reweighting which should contain systematic errors from Taylor expansion. Nevertheless, we see that the errors of the chemical potential from nonreweighted results overlap with those from the reweighted ones and the extrapolated crossing point of the two boundaries also agree with at albeit with larger errors.
We also compare our results to ones from a canonical ensemble study that uses staggered fermions [25] (Fig. 6 (right panel)). We find that the results are consistent with each other; however, our coexistence region is narrower. The discrepancy could come from the difference in the quark mass ( for our study compared to for the staggered study), and/or the fact that we use a different fermion formulation.
As a sanity check, we examine the seriousness of the potential sign problem. The average sign in Eq.(3.0) is plotted in Fig. 7 (left panel). We see that they range mostly from 0.6 to 0.2. Except for the point at and which has a 1.5 sigma away from zero, the others all have more than 3 sigmas above zero. In view of the fact that our results agree with those from the imaginary chemical potential study and the extrapolated boundaries meet at the expected at , we believe that the sign fluctuation in the reweighting is not a problem.
For the case, the phase diagram (see Fig. 3) is expected to be similar to the conjectured real QCD phase diagram. For nonzero quark mass, it also shows a crossover behavior at vanishing chemical potential, a critical point as well as a first order phase transition line. For this system, we scan three temperatures below : , and . The results are presented in Fig. 8: Given our limited statistics of 2000 configurations and the fact that the sign problem starts to set in below , we do not observe a clear signal for the Sshape structure within the scanned temperature and density range. This may be due to the fact that the transition is milder than the sensitivity of our data, or it may be due to the fact that the transition occurs at lower temperatures than the ones we studied – there is at least one study that claims that the critical point for occurs at temperatures below [35]. However, in view of the sign problem we encounter below , the claim needs to be scrutinized.
4.3 Phase diagram for
As we see from the above studies, the canonical approach at finite density is feasible and has been checked in the case. We would like to extend this study to the real world which contains two light quarks and one strange quark. However, simulating light quark masses requires larger volumes than the ones used in the present study. As a first attempt, we shall use the quark mass around physical strange quark mass as our input and apply the canonical ensemble method to the three degenerated flavor case which is expected to have similar features to the real QCD case. We are looking for the existence of a first order phase transition and its critical point at this stage, which may serve as an indication for the existence and approximate location of such a point in the real world. Following the methodology for two and four flavors, we fix the temperature at four different values and scan for the phase diagram by varying . The resulting baryon chemical potential is plotted in Fig. 9.
The presence of the Sshape structure at three temperatures below is a signal for a first order phase transition. Determination of the coexistence phase boundaries is carried out via the Maxwell construction. We plot the results from the Maxwell constructions in Fig. 9.
In Fig. 10 (left panel) we plot the boundary points as determined by the Maxwell construction. The phase boundaries intersect at the critical point which is located at and . In the right panel of Fig. 10, we plot the critical chemical potential as a function of temperature. On this graph the critical point is located at and .
Now that we have determined the phase boundaries in the temperature–density plot, we could map out the phase boundary in the temperature–chemical potential plot. The later is the usual phase diagram studied in the grand canonical ensemble. It is shown in Fig. 10 (right panel).
From our current simulations, we have found a first order phase transition and determined the location of the critical point for the three degenerate flavor case at intermediate quark mass. At the critical point, the critical temperature is determined as and the baryon chemical potential is . Given , we are supposed to see an Sshaped behavior at . However, we cannot observe a clear signal of it, it may be due to the fact that this temperature is too close to .
We show the average sign from reweighting (Eq. (3.0)) in Fig. 11 (left panel). Again, we see that most of them are more than 3 sigmas above zero. The only exceptions are those points at which are beyond the coexistence phase region. We also plot the average sign with no reweighting (Eq. (2.0)) in Fig. 11 (right panel). As expected, we find that they are in general larger in value than their respective ones with reweighting.
Before we conclude this section, we would like to discuss a potential contradiction between our results and those of de Forcrand and Philipsen [47, 48, 49]. These authors used 2+1 and 3 flavors staggered fermions and a Taylor expansion in to study the curvature of the critical surface at very light quark masses close to surface. After expanding to , they find that the critical surface bends such that the first order region shrinks at higher quark masses. This suggests an exotic scenario where there is no critical point at finite chemical potential [47, 48, 49]. This conclusion seems to contradict our results. However, since this conclusion is derived from studies at very low chemical potential, it is possible that the critical surface bends back at larger chemical potentials and the critical point reappears. In fact, an NJL model study seems to indicate that this scenario is actually realized [50].
5 Summary
In this study, we show that the canonical ensemble can be used to investigate the phase diagram at finite temperature and nonzero baryon density, at least on small lattices. In order to scan the QCD phase diagram at relatively large quark number (e.g. > 20), we developed the winding number expansion method to calculate the the projected determinant. It greatly expands our ability of scanning a broad region of the phase space and allows us to reach the coexistence phase.
We presented our results from simulations using two and four degenerate flavors of quarks. At finite volume, the first order phase transition at finite density appears as an Sshaped structure in the plot of chemical potential versus baryon number (baryon density). To determine the phase boundaries we use the Maxwell construction. For , the first order phase transition is observed at three different temperatures below the transition temperature at zero chemical potential. The phase boundary separating the hadron gas phase from the coexistence phase and the boundary between the plasma phase and the coexistence phase eventually cross at one point. For the four flavor case, this point should be the first order transition point at zero chemical potential. We have checked this numerically and found that the intersection of the boundaries indeed coincides with the the critical temperature at within errors. Furthermore, our results are consistent with those of the imaginary chemical potential approach in the grand canonical ensemble [45]. This is a powerful check for our method. For the two flavor case, we do not see any clear signal of a first order phase transition for temperature as low as .
The most interesting study is the three flavor case. We found signals of a first order phase transitions at three different temperatures and determined the phase boundaries for each temperature. The critical point located at and is obtained by an extrapolation.
We should point out that the present simulations are conducted on a small volumes, and with relatively large quark masses and a coarse lattice spacing ( fm). We plan to run simulations at lower quark masses using nHYP clover fermions.
Acknowledgments
This proceedings would not been possible without the work from my collaborators of the QCD collaboration. I am most grateful to Andrei Alexandru, KehFei Liu and Xiangfei Meng for their significant contributions and enjoyable collaboration. I would also like to thank P. de Forcrand, S. Ejiri, C. Gattringer, F. Karsch and M.P. Lombardo for their useful discussions and suggestions. This work is partially supported by DOE Grants DEFG0584ER40154.
References
 [1] F. Karsch and E. Laermann, (2003), arXiv:heplat/0305025.
 [2] S. Muroya, A. Nakamura, C. Nonaka, and T. Takaishi, Prog. Theor. Phys. 110, 615 (2003), arXiv:heplat/0306031.
 [3] K. Rajagopal and F. Wilczek, (2000), arXiv:hepph/0011333.
 [4] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Szabo, Nature 443, 675 (2006), arXiv:heplat/0611014.
 [5] M. Cheng et al., Phys. Rev. D74, 054507 (2006), arXiv:heplat/0608013.
 [6] HotQCD, C. E. Detar and R. Gupta, PoS LAT2007, 179 (2007), arXiv:0710.1655.
 [7] F. Karsch, PoS LAT2007, 015 (2007), arXiv:0711.0661.
 [8] RBC, F. Karsch, J. Phys. G35, 104096 (2008), arXiv:0804.4148.
 [9] Y. Aoki, Z. Fodor, S. D. Katz, and K. K. Szabo, Phys. Lett. B643, 46 (2006), arXiv:heplat/0609068.
 [10] Y. Aoki, S. Borsanyi, S. Durr, Z. Fodor, S. D. Katz, S. Krieg and K. K. Szabo, JHEP 0906, 088 (2009), [arXiv:0903.4155 [heplat]].
 [11] M. G. Alford, K. Rajagopal, and F. Wilczek, Nucl. Phys. B537, 443 (1999), arXiv:hepph/9804403.
 [12] F. R. Brown et al., Phys. Rev. Lett. 65, 2491 (1990).
 [13] M. Asakawa and K. Yazaki, Nucl. Phys. A504, 668 (1989).
 [14] A. Barducci, R. Casalbuoni, S. De Curtis, R. Gatto, and G. Pettini, Phys. Lett. B231, 463 (1989).
 [15] A. Barducci, R. Casalbuoni, S. De Curtis, R. Gatto, and G. Pettini, Phys. Rev. D41, 1610 (1990).
 [16] A. Barducci, R. Casalbuoni, G. Pettini, and R. Gatto, Phys. Rev. D49, 426 (1994).
 [17] J. Berges and K. Rajagopal, Nucl. Phys. B538, 215 (1999), arXiv:hepph/9804233.
 [18] A. M. Halasz, A. D. Jackson, R. E. Shrock, M. A. Stephanov, and J. J. M. Verbaarschot, Phys. Rev. D58, 096007 (1998), arXiv:hepph/9804290.
 [19] M. A. Stephanov, K. Rajagopal, and E. V. Shuryak, Phys. Rev. Lett. 81, 4816 (1998), arXiv:hepph/9806219.
 [20] M. A. Stephanov, K. Rajagopal, and E. V. Shuryak, Phys. Rev. D60, 114028 (1999), arXiv:hepph/9903292.
 [21] J. Engels, O. Kaczmarek, F. Karsch, and E. Laermann, Nucl. Phys. B558, 307 (1999), arXiv:heplat/9903030.
 [22] K. F. Liu, Chin. J. Phys. 38, 605 (2000).
 [23] K.F. Liu, Int. J. Mod. Phys. B16, 2017 (2002), arXiv:heplat/0202026.
 [24] V. Azcoiti, G. Di Carlo, A. Galante, and V. Laliena, JHEP 12, 010 (2004), arXiv:heplat/0409157.
 [25] P. de Forcrand and S. Kratochvila, Nucl. Phys. Proc. Suppl. 153, 62 (2006), arXiv:heplat/0602024.
 [26] K.F. Liu, (2003), arXiv:heplat/0312027.
 [27] A. Alexandru, M. Faber, I. Horvath, and K.F. Liu, Phys. Rev. D72, 114513 (2005), arXiv:heplat/0507020.
 [28] A. Alexandru, A. Li, and K.F. Liu, PoS LAT2007, 167 (2007), arXiv:0711.2678.
 [29] X.f. Meng, A. Li, A. Alexandru, and K.F. Liu, PoS LATTICE2008, 032 (2008), arXiv:0811.2112.
 [30] J. Danzer, C. Gattringer, and L. Liptak, PoS LAT2009, 185 (2009), arXiv:0910.3541.
 [31] C. Gattringer and L. Liptak, (2009), arXiv:0906.1088.
 [32] J. Danzer and C. Gattringer, Phys. Rev. D78, 114506 (2008), arXiv:0809.2736.
 [33] A. Li, A. Alexandru, and K.F. Liu, PoS LAT2006, 030 (2006), arXiv:heplat/0612011.
 [34] A. Li, A. Alexandru, and K.F. Liu, PoS LAT2007, 203 (2007), arXiv:0711.2692.
 [35] S. Ejiri, Phys. Rev. D78, 074507 (2008), arXiv:0804.3227.
 [36] Y. Iwasaki, K. Kanaya, S. Kaya, and T. Yoshie, Phys. Rev. Lett. 78, 179 (1997), arXiv:heplat/9609022.
 [37] CPPACS, A. Ali Khan et al., Phys. Rev. D63, 034502 (2001), arXiv:heplat/0008011.
 [38] CPPACS, A. Ali Khan et al., Phys. Rev. D64, 074510 (2001), arXiv:heplat/0103028.
 [39] S. A. Gottlieb, W. Liu, D. Toussaint, R. L. Renken, and R. L. Sugar, Phys. Rev. D35, 2531 (1987).
 [40] S. Duane, A. D. Kennedy, B. J. Pendleton, and D. Roweth, Phys. Lett. B195, 216 (1987).
 [41] A. D. Kennedy, I. Horvath, and S. Sint, Nucl. Phys. Proc. Suppl. 73, 834 (1999), arXiv:heplat/9809092.
 [42] M. A. Clark and A. D. Kennedy, Phys. Rev. Lett. 98, 051601 (2007), arXiv:heplat/0608015.
 [43] B. Joo, I. Horvath, and K. F. Liu, Phys. Rev. D67, 074505 (2003), arXiv:heplat/0112033.
 [44] A. Li, A. Alexandru, and K.F. Liu, (in preparation) .
 [45] M. D’Elia and M.P. Lombardo, Phys. Rev. D67, 014505 (2003), arXiv:heplat/0209146.
 [46] O. Philipsen, Prog. Theor. Phys. Suppl. 174, 206 (2008), arXiv:0808.0672.
 [47] P. de Forcrand and O. Philipsen, PoS LATTICE2008, 208 (2008), arXiv:0811.3858.
 [48] P. de Forcrand and O. Philipsen, JHEP 11, 012 (2008), arXiv:0808.1096.
 [49] P. de Forcrand and O. Philipsen, PoS LAT2006, 130 (2006), arXiv:heplat/0611027.
 [50] J.W. Chen, K. Fukushima, H. Kohyama, K. Ohnishi, and U. Raha, Phys. Rev. D80, 054012 (2009), arXiv:0901.2407.