First order transition regions in the quark masses and chemical potential parameter space of QCD
We investigate the phase transitions of -flavor QCD, where two light flavors and massive flavors exist, aiming to understand the phase structure of (2+1)-flavor QCD. Performing simulations of 2-flavor QCD with improved staggered and Wilson fermions and using the reweighting method, we calculate probability distribution functions in the many-flavor QCD. Through the shape of distribution functions, we determine the critical surface terminating first order phase transitions in the parameter space of the light quark mass, heavy quark mass and the chemical potential, and find that the first order region becomes larger with . We then study the critical surface at finite density for large and the first order region is found to become wider with the increasing chemical potential. On the other hand, the light quark mass dependence of the critical mass of heavy quarks seems weak in the region we investigated. The result of this weak dependence suggests that the critical mass of heavy quark remains finite in the chiral limit of 2-flavors and there exists a second order transition region on the line of the 2-flavor massless limit above the tri-critical point. Moreover, we extend the study of 2-flavor QCD at finite density to the case of a complex chemical potential and investigate the singularities where the partition function vanishes, so-called Lee-Yang zeros. The plaquette effective potential is computed in the complex plane. We find that the shape of the effective potential changes from single-well on the real axis to double-well at large imaginary chemical potential and the double-well potential causes the singularities.
First order transition regions in the quark masses and chemical potential parameter space of QCD
Department of Physics, Saga University, Saga 840-8502, Japan
QCD at high temperature and density has rich phase structure, and the nature of the phase transition changes depending on the quark mass and the number of flavors. The quark mass dependence of the QCD phase transition is important not only for the theoretical interest in the phase structure of QCD with various quark masses but also for the investigation of finite density QCD or QCD near the chiral limit, in which the numerical studies are difficult. From the numerical study of (2+1)-flavor QCD, the QCD phase transition is considered to be crossover at low density but is expected to turn into a first order transition at high density. Finding the endpoint of first order transition is then crucial for establishing the above expectation. Unfortunately, it is still extremely difficult to simulate quarks with physical mass at high density. But accessing the endpoint becomes easy if one extends QCD in appropriate directions, and importantly if such an extension is smooth, one can extract information on original QCD by suitable extrapolations [1, 2].
Moreover, the nature of the phase transition in the chiral limit of 2-flavor QCD is a long-standing problem. The standard expectation is of second order, but the order is not conclusive due to the difficulty of the numerical study in the chiral limit. One of the good approaches is to investigate the boundary of the first order transition region. If the critical value of the strange quark mass does not go to infinity in the up and down quarks massless limit, the transition of the massless 2-flavor QCD is not of first order. However, it is difficult to study it, because the first order region is very small in (2+1)-flavor QCD and simulations with very small mass are required.
In Ref. , the boundary of the first order region is studied in many-flavor QCD, motivated by a feasibility study of the electroweak baryogenesis in technicolor theories constructed by many-flavor gauge theory. They studied QCD with two degenerate light quarks of the mass and the chemical potential and massive quarks with and , and found by measuring probability distribution function that the critical massive quark mass becomes larger as increases. Therefore, the investigation of the critical line becomes easier as increases. This extension of QCD is also useful for the study of the phase transition of massless 2-flavor QCD. If the critical mass for the massive flavors remains finite in the massless limit of two light flavors, it gives a strong support for the second order transition. A similar approach has been tried in Ref. .
This paper consists of two parts. The first part deals with the light quark mass dependence of the critical massive quark mass separating the first order and crossover regions in -flavor QCD. We perform simulations with 2 flavors of improved Wilson quarks. The effect from the dynamical -flavors are added by a reweighting method assuming that the -flavors are heavy. We then investigate the critical heavy quark mass as a function of the light quark mass through the shape of distribution functions and discuss the nature of phase transitions in the chiral limit of 2-flavor QCD.
The second part is the chemical potential dependence. In particular, we extend the real chemical potential to complex value. The phase transition of 2-flavor QCD with finite mass at is crossover. But, when the chemical potential is introduced, the shape of the distribution changes and the distribution function will become first-order-transition-like. For the complex , this property make a singularity, which is called “Lee-Yang zero” . Lee and Yang proposed the method to investigate the nature of phase transitions from the singularities in the complex parameter plane, and applications to finite density QCD have been discussed in Refs. [6, 7].
In the next section, we explain our method to identify the nature of phase transitions via the distribution function. We then argue the light quark mass dependence of the endpoint of the first order transition in Sec. 3. The singularities in the complex plane are discussed in Sec. 4. Conclusions are given in Sec. 5.
2 Critical point by a histogram method
We study QCD with two degenerate light quarks of the mass , the chemical potential and heavy quarks. We define the probability distribution function of average plaquette value,
where and are the actions of gauge and quark fields, respectively, and is the quark matrix. For simplification, we adopt which does not depend on explicitly. is the number of sites and is the simulation parameter. defined by is Wilson loop for the standard plaquette gauge action. is the delta function, which constrains the operator to be the value of . We moreover define the effective potential, .
Denoting the potential of 2-flavor at by , that of -flavor is written as
where and means the ensemble average over 2-flavor configurations generated at , and vanishing . is the simulation point, which may differ from in this method.
At a first order transition point, shows a double-well shape as a function of , and equivalently the slope of the potential shows an S-shape. Since appears only in the linear term of on the right hand side of Eq. (2.0), the shape of the slope is independent of , i.e. a change of just shifts the overall constant of the slope of . Although must be adjusted to the first order transition point to observe the double-well potential, the fine tuning is not necessary if we investigate the slope.
The derivative can be measured easily from the peak position of the plaquette histogram . When one performs a simulation at , the slope is vanishing at the minimum of , and the value of at the minimum can be estimated by approximately. Hence, we obtain at by
We therefore focus on the slope of the effective potential to identify the nature of transitions.
3 Light quark mass dependence of the critical heavy quark mass
In this section, we discuss the light quark mass dependence of the critical hopping parameter of heavy quarks terminating the first order transition. It is important to investigate whether the critical hopping parameter vanishes at a finite light quark mass or not. If the critical line crosses the line of , the 2-flavor chiral limit is in the first order region, where is the hopping parameter of heavy-flavors being proportional to the inverse mass. Restricting the calculation to the heavy quark region near , the second determinant for flavors in Eq. (2.0) is approximated by the leading order of the hopping parameter expansion,
for the standard Wilson quark action . and are the real and imaginary part of the Polyakov loop, respectively. For improved gauge actions such as , additional terms must be considered, where is the improvement coefficient and . However, since the improvement term does not affect the physics, we will cancel these terms by a shift of the coefficient . It is shown in Ref.  that the hopping parameter expansion is applicable for large .
Denoting , we obtain for with
where is given by the Polyakov loop and is independent of . The plaquette term does not contribute to and can be absorbed by shifting for Wilson quarks. Moreover, one can deal with the case with non-degenerate masses by adopting for the Wilson quark action or for the staggered quark action. Thus, the choice of the quark action is not important. In the following, we discuss the mass dependence of through the parameter .
We perform simulations of 2-flavor QCD with the clover-improved Wilson quark and Iwasaki gauge actions111The improvement parameters are almost the same as Ref. . We take four different values of light quark masses ranging from to . The corresponding ratio of pseudo-scalar and vector meson masses is and for and , respectively. The data are taken at about 30 values of around the pseudo-critical point at each , and at each simulation point 10,000 to 40,000 trajectories have been accumulated.
In the calculation of , we use the delta function approximated by , where is selected consulting the resolution and the statistical error. Because is independent of , we mix all data obtained at different as was done in Ref. . The histograms of plaquette value are plotted in the left panel of Fig. 1. The results for at are shown in the right panel of Fig. 1 for – at interval of . A rapid increase is observed around , and the gradient becomes larger as increases.
The derivative is calculated by fitting to a -order polynomial of in an appropriate fit range. Using the equation,
we compute for each . The results for various are shown in Fig. 2, where is computed by Eq. (2.0). The shape of is independent of because is -independent. The first derivative is the monotonically increasing function of when is small, indicating that the transition is crossover. However, the shape of turns into an S-shaped function at , which means that the system undergoes first order transition. The same analysis has been done in Ref.  using the p4-improved staggered fermion action for 2-flavor QCD with . The result of the critical value of at which the first order transition appears, , is about 0.06. The difference may be caused by the lattice discretization error due to small . We have defined the parameter for the Wilson quark. Then, the critical point corresponding decreases as with , and the truncation error from the higher order terms of the hopping parameter expansion in becomes smaller as increases.
The remarkable point of this study is light quark mass dependence. In Fig. 3, we plot the results of the critical value as functions of (left) and the quark mass defined by the PCAC relation (right). The slope of is computed fitting the data with or order polynomials. The open symbols are the results by order, and the filled symbols are those by order. The difference is taken as the systematic error. The red lines are the upper and lower bounds of the critical for the system with massive quarks and no light quarks, which is determined by observing the hysteresis curves of the plaquette and Polyakov loop. It is found from these figures that the light quark mass dependence of the critical mass of heavy quarks is very small in the region we computed. The weak dependence suggests that the critical mass of heavy quark remains finite in the chiral limit of 2-flavors. Thus, the sign of the first order transition in 2-flavor QCD is not shown.
4 Singularities in the complex chemical potential plane
Let us turn to finite density QCD. The distribution function of 2-flavor QCD at finite density and the appearance of the double-well potential are discussed in Ref. , and the boundary of the first order transition region of -flavor QCD at finite density is computed in Ref.  using the 2-flavor QCD configurations generated with the p4-improved staggered quark action in Ref. . The first order region is found to become larger as increases.
In this section, we extend the potential analysis to the complex chemical potential, . The numerical study of the singularities of has a potential danger for large when we use the reweighting method . We estimate the position of from the probability distribution of the complex phase, since vanishes when the distribution function has two peaks or more and the contributions from these peaks cancel each other . In this study, we investigate the plaquette distribution function and complex phase average with constraining the plaquette value. To avoid the sign problem, the Gaussian approximation  is applied 222Preliminary results are presented in Ref. .
We compute the probability distribution function of the plaquette ,
We denote . The normalized partition function is rewritten as
Here, is the reweighting factor for finite defined by Eq. (2.0) with and for -flavors. This is independent of and can be measured at any . In this study, all simulations are performed at and the effect of finite is introduced though the operator measured on the configurations generated by the simulations at .
Because QCD has time-reflection symmetry, the partition function is invariant under a change from to , i.e. . Moreover, the quark determinant satisfies . From these equations, we get
This indicates that is real valued function in the case of real , i.e. . Then, the plaquette probability distribution is also real valued. However, once the imaginary part of becomes nonzero, is not a real number any more. We thus write the partition function,
This complex phase vanishes at , and is monotonic function of at small if we define the complex phase by a Taylor expansion of .
If becomes a double-peaked function which has two peaks of equal height at and , two phases coexist like first order phase transitions. When changing , the expectation value of plaquette changes discontinuously form to at that point. And then, the partition function can be approximated by . When the difference between and is equal to with an integer , the partition function will be vanishing, i.e. Lee-Yang zeros appear. The first-order-transition-like point in the complex plane corresponds to “Stokes line” in the infinite volume limit. We define the effective potential as and investigate the position of this would-be Stokes line in the complex plane.
We calculate these three quantities, , and by Monte-Carlo simulations. Because the exact calculation of the quark determinant is difficult except on small lattices, we estimate the quark determinant from the data of Taylor expansion coefficients of up to around obtained by a simulation of 2-flavor QCD with p4-improved staggered quarks in Ref. . The truncation error has been discussed in Ref. . Notice that the complex phase is not restricted to the range from to because we define by the Taylor expansion of .
To compute at finite , we discuss the sign problem. We denote the quark determinant as . Histograms of seem to be well-approximated by Gaussian functions. (See Fig. 1 (left) in Ref. .) Here, we perform a cumulant expansion,
For the case that the distribution of is of Gaussian, the terms vanish for in this equation. Moreover, since and , this expansion can be regarded as a power expansion in . If the convergence of the expansion Eq. (4.0) is good, we can extract the phase factor from easily and the sign problem in is eliminated. We deal with the first two terms, i.e. and , assuming the Gaussian distribution. The correlation terms between and are also neglected as a first step. Because we calculate the expectation value with fixed and the values of and are strongly correlated, the may be small once is fixed. Then, .
We discuss instead of itself, again. at different can be estimated by the equation,
under the parameter change from to .
If the effective potential is a double-well function of having two minima at and and one maximum at the middle , is an S-shaped function and vanishes three times at , and . The condition,
is satisfied when at the transition point. In such a case, there exists a region of where the derivative of is negative. The results of the first derivative of are shown in Fig. 4 for various with (left) and (right), where . We measured them at the peaks of the plaquette histograms and interpolated the data by a cubic spline method. In the region of large , becomes an S-shaped function. The result of is plotted in Fig. 5 (left) for . It is a monotonically increasing function of .
We then investigate the boundary at which changes to an S-shaped function from a monotonic function. The boundary is shown as a red line in Fig. 5 (right). Above this line, the region of where appears. The constant part of is changed by varying . We then adjust such that the depths of two minima of are the same using Eq. (4.0). If the dashed lines in Fig. 4 move to the horizontal axis by changing , Eq. (4.0) is satisfied.
Contour plots of the at which are shown by blue curves in Fig. 5 (right). The values of are indicated near the blue lines. Along the line for each , the double-well potential appears, and the contour curve is expected to turn into the Stokes line in the infinite volume limit. For finite volume, at points on that line, where the phase cancelation occurs, Lee-Yang zeros appear. This result indicates the existence of singularities in the region of large as well as the region of large , and the boundary in the complex plane is closer to the origin than that (the square symbol) on the real axis.
We studied the distribution function and effective potential of QCD with two light flavors and massive flavors, aiming to understand phase structure of 2 and (2+1)-flavor QCD. Through the shape of the distribution function, we investigated the critical surface separating the first order transition and crossover regions in the parameter space of the light quark mass, heavy quark mass and the chemical potential. It is found that the critical mass becomes larger with and the first order region becomes wider with the increasing chemical potential for large . If (2+1)-flavor QCD has the same property, this gives the strong evidence for the existence of the critical point at finite density in the real world.
The nature of the chiral phase transition in the 2-flavor massless limit is still open question. To study the chiral limit, we investigated the light quark mass dependence. If the transition is of first order, the critical vanishes before going to the 2-flavor massless limit. But, the critical does not show such a behavior in the region we investigated. This implies that the critical mass of heavy quark remains finite even in the chiral limit of 2-flavors and 2-flavor QCD near the chiral limit is not in the first order transition region.
We then discussed 2-flavor QCD with the complex by a numerical simulation. We found that there is a region where the plaquette distribution function has two peaks, suggesting the existence of singularities characterized by , in the region of large as well as of large . Combined with the analytical study of the complex chemical potential, the distribution of the singularities may provide important information about the QCD phase transition.
-  Ch. Schmidt, et al., Nucl. Phys. B (Proc. Suppl.) 119, 517 (2003); F. Karsch, et al., Nucl. Phys. B (Proc. Suppl.) 129, 614 (2004); S. Ejiri, et al., Prog. Theor. Phys. Suppl. 153, 118 (2004).
-  P. de Forcrand and O. Philipsen, Nucl. Phys. B 673, 170 (2003); JHEP 0701 077 (2007).
-  S. Ejiri and N. Yamada, Phys. Rev. Lett. 110, 172001 (2013).
-  C. Bonati, P. de Forcrand, M. D’Elia, O. Philipsen, and F. Sanfilippo, Phys. Rev. D 90, 074030 (2014).
-  C.N. Yang and T.D. Lee, Phys. Rev. 87, 404 (1952); T.D. Lee and C.N. Yang, Phys. Rev. 87, 410 (1952); M.E. Fisher, Phys. Rev. Lett. 40, 1610 (1978).
-  M.A. Stephanov, Phys. Rev. D 73, 094508 (2006).
-  S. Ejiri, Y. Shinno and H. Yoneyama, Prog. Theor. Exp. Phys. 2014, 083B02 (2014).
-  S. Ejiri, Phys. Rev. D 77, 014508 (2008).
-  S. Ejiri and H. Yoneyama, PoS(LATTICE 2009) 173 (2009).
-  H. Saito, S. Ejiri, S. Aoki, K. Kanaya, Y. Nakagawa, H. Ohno, K. Okuno, T. Umeda (WHOT-QCD Collaboration), Phys. Rev. D 89, 014508 (2014).
-  S. Ejiri, Y. Maezawa, N. Ukita, S. Aoki, T. Hatsuda, N. Ishii, K. Kanaya, and T. Umeda (WHOT-QCD Collaboration), Phys. Rev. D 82, 014508 (2010).
-  S. Ejiri, Phys. Rev. D 73, 054502 (2006).
-  C.R. Allton, M. Döring, S. Ejiri, S.J. Hands, O. Kaczmarek, F. Karsch, E. Laermann and K. Redlich, Phys. Rev. D 71, 054508 (2005).