# Critical point of QCD from lattice simulations in the canonical ensemble

###### Abstract

A canonical ensemble algorithm is employed to study the phase diagram of QCD using lattice simulations. We lock in the desired quark number sector using an exact Fourier transform of the fermion determinant. We scan the phase space below and look for an S-shape structure in the chemical potential, which signals the coexistence phase of a first order phase transition in finite volume. Applying Maxwell construction, we determine the boundaries of the coexistence phase at three temperatures and extrapolate them to locate the critical point. Using an improved gauge action and improved Wilson fermions on lattices with a spatial extent of and quark masses close to that of the strange, we find the critical point at and baryon chemical potential .

###### pacs:

11.15.Ha, 11.30.Rd^{†}

^{†}preprint: UK/11-02{CJK*}

UTF8 \CJKfamilygbsn

QCD Collaboration

QCD is expected to have a rich phase diagram at finite temperature and finite density. Current lattice calculations have shown that the transition from the hadronic phase to QGP phase is a rapid crossover Aoki:2009sc (); Karsch:2008fe (). For large baryon chemical potential and very low temperature, a number of models suggest that the transition is a first order. If this is the case, when the chemical potential is lowered and temperature raised, this first order phase transition is expected to end as a second order phase transition point — the critical point. However, lattice QCD simulations with chemical potential are difficult due to the notorious “sign problem”. The majority of current simulations are focusing on small chemical potential region where the “sign problem” appears to be under control. Up to now, all the or 2+1 simulations are based on the grand canonical ensemble (, as parameters) with staggered fermions. The results from the multi-parameter reweighting Fodor:2004nz (), Taylor expansion with small Allton:2005gk (); Gavai:2008zr () and the curvature of the critical surface deForcrand:2008vr () are not settled and need to be cross-checked. Even the existence of the critical point is in question deForcrand:2008vr (). We employ an algorithm, which is not restricted to small chemical potential because of the mitigation of the sign problem under the current parameter settings, to study this problem.

In this letter, we adopt an exact Monte Carlo algorithm Liu:2003wy (); Alexandru:2005ix (); Alexandru:2007bb () based on the canonical partition function Engels:1999tz (); Liu:2000dj (); Liu:2002qr (); Azcoiti:2004ri (); deForcrand:2006ec (); Kratochvila:2005mk () which is designed to alleviate the determinant fluctuation problems. As it turns out, the sign fluctuations are not serious on the lattices used in the present study, as we shall see later. In the canonical ensemble simulations in finite volume, the coexistence phase of a first order phase transition has a characteristic S-shape as a function of density due to the surface tension. This finite-volume property has been exploited successfully to identify the phase boundaries via the Maxwell construction in studies of phase transition with the staggered fermions deForcrand:2006ec (); Kratochvila:2005mk () and clover fermions Li:2010qf () for the case which is known to have a first order phase transition at . In these benchmark studies the boundaries were identified at three temperatures below , and they were extrapolated in density and temperature to show that the intersecting point indeed coincides with the independently identified first order transition point at and Li:2010qf (). In view of the success of the study, we extend this method to the more realistic case Li:2009ju (); Li:2010dya (). Although the real world contains two light quarks and one heavier strange quark, the three degenerate flavor case has a similar phase structure. Our primary goal in this study is to determine whether a first order phase transition exists for and where the critical point is located.

With the aid of recently developed matrix reduction technique Alexandru:2010yb (); Nagata:2010xi (); Gattringer:2009wi (), we scan the chemical potential as a function of baryon number for four temperatures below which is determined at zero chemical potential, and we observe clear signals for a first order phase transition for temperatures below . The phase boundaries of the coexistence phase are determined and then extrapolated in temperature and density to locate the critical point at and . Our results are based on simulations on lattices with clover fermion action with quark masses which correspond to the pion mass from for the lowest temperature to for the highest temperature.

The canonical partition function in lattice QCD can be derived from the fugacity expansion of the grand canonical partition function,

(1) |

where is the net number of quarks (number of quarks minus the number of anti-quarks) and is the canonical partition function. Using the fugacity expansion, it can be shown that the canonical partition function can be written as a Fourier transform of the grand canonical partition function,

(2) |

upon introducing an imaginary chemical potential . After integrating out the fermionic part in Eq. (2), we get an expression

(3) |

where

(4) |

is the projected determinant for the fixed net quark number . is the number of flavors. We shall use the recently developed matrix reduction technique to compute the projected determinant exactly Alexandru:2010yb ().

Using charge conjugation symmetry, one can show that is real, but not necessarily positive. Due to the sign fluctuation, there can potentially be a sign problem at large quark number and low temperature. For more detailed discussion about the properties of the canonical ensemble, we refer the reader to Ref. Li:2010qf (). To simulate Eq. (3) dynamically, we rewrite canonical partition function as

(5) |

where

(6) |

Our strategy to generate an ensemble is to employ Metropolis accept/reject method based on the weight and fold the phase factor into the measurements. In short, during the simulation, the candidate configuration is “proposed” by the standard Hybrid Monte Carol algorithm and then an accept/reject step is used for the correct probability. Note the two-step simulation with HMC and accept/reject based on reduces the fluctuation problem Alexandru:2005ix () and accept/reject step based on the exact projected determinant ensures that the simulation remains in the specific canonical sector with quark number .

The lattice spacing and the pion mass are determined by using dynamically generated ensembles on lattices for each . To locate the pseudo critical temperature , we varied to look for the peak of the Polyakov loop susceptibility. We run simulations for five different volumes () and found that the peak of the susceptibility hardly depends on the volume. This is consistent with the finding on large volumes and physical quark masses that the finite temperature transition for the case is a crossover at zero chemical potential Aoki:2009sc (); Karsch:2008fe ().

To determine the location of the phase transition at non-zero baryon density, we pick four temperatures below (, , , ) and vary the net quark number from 3 to 54 in steps of 3 (for fractional baryon number the partition function vanishes). This corresponds to the baryon number from 1 to 18 and a density between that of the nuclear matter and 18 times of that. The chemical potential is calculated and plotted as a function of the net baryon number . In the canonical ensemble, the baryon chemical potential is calculated by taking the difference of the free energy after adding one baryon, i.e.

(7) |

where

(8) |

is measured in the ensemble with baryon number and in Eq. (7) stands for the average over the ensemble generated with the measure .

As a first check, we examine the magnitude of the sign fluctuations. The average sign in Eq. (6) appears in the denominator of Eq. (7) and can lead to a sign problem when its error bars overlap with zero. This quantity is plotted in Fig. 2 for the highest and lowest temperatures. We see that all of them are more than above zero. This result is better than the previous ones based on the winding number expansion method Li:2010qf (); Li:2010dya (), presumably due to the adoption of the exact projection of the determinant Alexandru:2010yb (). Thus, we believe that the sign fluctuations are not a problem for this study.

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 expected canonical ensemble phase diagram in Fig. 4 in contrast to that in the grand canonical ensemble in Fig. 1. The first order phase transition line in the grand canonical diagram becomes a phase coexistence region in the - diagram of the canonical ensemble, which has two boundaries that separate it from the pure phases. The two boundaries will eventually meet at one point. This point is the critical point at nonzero baryon chemical potential.

Once one enters the coexistence region in a finite volume, the contribution from the surface tension causes the appearance of a “double-well” in the effective free energy whose derivative with respect to density leads to an S-shaped behavior in the chemical potential versus baryon number plot Ejiri:2008xt (). However, in the thermodynamic limit, the surface tension contribution goes away since it is a surface term while the free energy scales with the volume; the chemical potential will then stay constant in the coexistence phase region. The behavior of the baryon chemical potential in the thermodynamic limit is shown as an inset in Fig. 4. and mark the lower and upper boundaries of the coexistence phase at a given temperature below .

Our results for the baryon chemical potential are presented in Fig. 5 for four different temperatures below . Statistical errors are estimated from the jackknife method. It is clear that the chemical potential exhibits an “S-shaped” wiggle for between and . To identify the boundaries of the mixed-phase region and the coexistence baryon chemical potential, we rely on the Maxwell construction: the coexistence chemical potential is the one that produces equal areas between the curve of the chemical potential as a function of and the constant line which intersects with at and . This procedure was used in studies with staggered fermions deForcrand:2006ec (); Kratochvila:2005mk () and Wilson-clover fermions Li:2010qf () in this context for the case.

We carried out the Maxwell constructions for the three temperatures at , and . We could not do it for the 0.93 case, as the wiggle there, if present, is not statistically significant. The results are presented in Fig. 3. Having determined and for three temperatures, we plot the boundaries of the coexistence region and perform an extrapolation in and to locate the intersection of the two boundaries. To determine the crossing point, we perform a simultaneous fit of the boundary lines using a even polynomial in baryon density. We use an even polynomial since is an even function of . The phase boundaries and their extrapolations are plotted in Fig. 6. We find the intersection point at and .

Using the coexistence chemical potential, one can map out the phase diagram in the grand canonical ensemble as shown in Fig. 7. Note that, the region of coexistence phase becomes a curved transition line separating two the phases as we expected. In this way, we locate the critical point in the grand canonical ensemble at critical temperature and baryon chemical potential . Using the lattice spacing in our simulation, we convert its location in physical units to be and .

In conclusion, we have applied a canonical ensemble algorithm previously tested on the to the more relevant case and located the first order phase transition as signaled by the S-shape structure in the plane for several temperatures below . The Maxwell construction was employed to identify the boundaries of the coexistence phase and we extrapolated them to locate the critical point at and . We should point out that the present work is carried out on a relatively small volume with spatial extent of and for three degenerate quark flavors with their masses similar to that of the strange quark. Quark mass for this system acts like the magnetic field for spin systems which weakens the phase transition. Since the finite temperature transition is first order for massless quarks Pisarski:1983ms () and the present critical point is at a relatively large for quark masses around the strange, one expects that the critical point for the more realistic flavor case with light quarks to be somewhere in between. This expectation is based on the assumption that there is a critical surface which grows continuously from the critical line at into finite . The critical line is the one that separates the first order phase transition region at some finite temperature with small quark masses and the crossover region with intermediate masses (including the physical ones) on the plane of the Columbia plot. This assumption is challenged by recent studies of the critical surface near the critical line (See Fig.1 in references deForcrand:2008vr (); deForcrand:2008zi ()) which suggest that the first order region shrinks with increasing chemical potential and, therefore, there might not be a critical point for physical quark masses. To address this issue, future simulations will study the quark mass dependence of the critical point and the existence of the critical point needs to be checked on lattices with higher cutoffs and larger volumes.

This work is partially supported by DOE grants DE-FG05-84ER40154, DE-FG02-05ER41368 and DE-FG02-95ER-40907. We wish to thank P. de Forcrand, A. Kennedy and S. Chandrasekharan for useful discussions.

## References

- (1) Y. Aoki et al., JHEP 06, 088 (2009), arXiv:0903.4155 [hep-lat]
- (2) F. Karsch (RBC), J. Phys. G35, 104096 (2008), arXiv:0804.4148 [hep-lat]
- (3) Z. Fodor and S. Katz, JHEP 0404, 050 (2004), arXiv:hep-lat/0402006 [hep-lat]
- (4) C. Allton et al., Phys.Rev. D71, 054508 (2005), arXiv:hep-lat/0501030 [hep-lat]
- (5) R. Gavai and S. Gupta, Phys.Rev. D78, 114503 (2008), arXiv:0806.2233 [hep-lat]
- (6) P. de Forcrand and O. Philipsen, JHEP 11, 012 (2008), arXiv:0808.1096 [hep-lat]
- (7) K.-F. Liu, Edinburgh 2003, QCD and numerical analysis III (Springer-Verlag), 101(2005), arXiv:hep-lat/0312027
- (8) A. Alexandru, M. Faber, I. Horvath, and K.-F. Liu, Phys. Rev. D72, 114513 (2005), arXiv:hep-lat/0507020
- (9) A. Alexandru, A. Li, and K.-F. Liu, PoS LAT2007, 167 (2007), arXiv:0711.2678 [hep-lat]
- (10) J. Engels et al., Nucl. Phys. B558, 307 (1999), arXiv:hep-lat/9903030
- (11) K. F. Liu, Chin. J. Phys. 38, 605 (2000)
- (12) K.-F. Liu, Int. J. Mod. Phys. B16, 2017 (2002), arXiv:hep-lat/0202026
- (13) V. Azcoiti et al., JHEP 12, 010 (2004), arXiv:hep-lat/0409157
- (14) P. de Forcrand and S. Kratochvila, Nucl. Phys. Proc. Suppl. 153, 62 (2006), arXiv:hep-lat/0602024
- (15) S. Kratochvila and P. de Forcrand, PoS LAT2005, 167 (2006), arXiv:hep-lat/0509143
- (16) A. Li, A. Alexandru, K.-F. Liu, and X. Meng, Phys. Rev. D82, 054502 (2010), arXiv:1005.4158 [hep-lat]
- (17) A. Li, A. Alexandru, X. Meng, and K.-F. Liu (chi QCD), Nucl. Phys. A830, 633c (2009), arXiv:0908.1155 [hep-lat]
- (18) A. Li, PoS LAT2009, 011 (2009), arXiv:1002.4459 [hep-lat]
- (19) A. Alexandru and U. Wenger, Phys. Rev. D83, 034502 (2011), arXiv:1009.2197 [hep-lat]
- (20) K. Nagata and A. Nakamura, Phys. Rev. D82, 094027 (2010), arXiv:1009.2149 [hep-lat]
- (21) C. Gattringer and L. Liptak, Phys. Lett. B697, 85 (2011), arXiv:0906.1088 [hep-lat]
- (22) S. Ejiri, Phys. Rev. D78, 074507 (2008), arXiv:0804.3227 [hep-lat]
- (23) R. D. Pisarski and F. Wilczek, Phys. Rev. D29, 338 (1984)
- (24) P. de Forcrand and O. Philipsen, PoS LATTICE2008, 208 (2008), arXiv:0811.3858 [hep-lat]