Ab initio study of the thermodynamics of quantum chromodynamics on the lattice at zero and finite densities

Ab initio study of the thermodynamics of quantum chromodynamics on the lattice at zero and finite densities Shinji Ejiri,1,∗ Kazuyuki Kanaya,2,∗ and Takashi Umeda3,∗ for the WHOT-QCD Collaboration§ 1Graduate School of Science and Technology, Niigata University, Niigata 950-2181, Japan 2 Faculty of Pure and Applied Sciences, University of Tsukuba, Tsukuba, Ibaraki 305-8571, Japan 3 Graduate School of Education, Hiroshima University, Higashihiroshima, Hiroshima 739-8524, Japan §Sinya Aoki, Shinji Ejiri, Tetsuo Hatsuda, Noriyoshi Ishii, Kazuyuki Kanaya, Yoshiyuki Nakagawa, Yu Maezawa, Hiroshi Ohno, Hana Saito, Naoya Ukita, and Takashi Umeda. ∗E-mail: ejiri@muse.sc.niigata-u.ac.jp; kanaya@ccs.tsukuba.ac.jp; tumeda@hiroshima-u.ac.jp


Introduction
Temperature, T , and density are controllable parameters in a system. At sufficiently high T , we expect that confinement is violated and chiral symmetry is recovered because the effective coupling at the thermal energy scale becomes small due to the asymptotic freedom of quantum chromodynamics (QCD). At this point, systems with quarks and gluons will form a colored plasma state called "quark-gluon plasma" (QGP). Although humankind has never experienced QGP, it is expected to play an important role in the creation of matter during the early development of the Universe. Furthermore, QGP is considered to be observed by relativistic heavy-ion collision experiments at RHIC and LHC [1].
As in the high-temperature case, we expect deconfinement at sufficiently high density, because the average distance between quarks becomes small at this point and the properties of the system will be dominated by the asymptotic freedom. The density is controlled by the chemical potential μ. At very large μ and low T , we expect a BCS-like state called a "color superconductor" due to the attractive interaction between quarks. At lower densities around the nuclear density, we expect a nuclear fluid state, which may appear around the core of neutron stars. We thus expect a rich phase structure in PTEP 2012, 01A104 S. Ejiri et al.  = m u = m d and the strange quark mass m s . The top-right corner corresponds to the quenched limit of QCD. Lattice simulations with improved staggered quarks suggest that the physical point is located in the crossover region. The phase diagram is not fully established yet. See Refs. [3][4][5] for discussions and caveats.
QCD as a function of T and μ. See, e.g., Ref. [2] for a recent review. A prospective phase diagram is shown in Fig. 1.
When we vary the quark masses off the physical point, the nature of the quark matter may be different as a result. The usual expectation for the order of the finite-temperature QCD transition at μ = 0, based on effective model studies and lattice simulations, is summarized in Fig. 2. The details of the phase diagram as well as the nature of the quark matter at finite T and μ are, however, not well clarified yet. Because the issue is essentially nonperturbative, numerical studies on the lattice are the only systematic way to investigate the phase structure directly from the first principles of QCD [3][4][5][6][7][8][9].
Most lattice studies of hot/dense QCD have been done with computationally less demanding staggered-type lattice quarks [3][4][5][6][7][8][9]. In particular, in the study of the equation of state (EOS), extrapolations to the physical point and to the continuum limit have been achieved only with staggered-type quarks. However, theoretical bases such as locality and universality are not well established with them [3][4][5]. Therefore, to evaluate the effects of lattice artifacts and thus to obtain reliable 2/35 PTEP 2012, 01A104 S. Ejiri et al. predictions to be compared with experiment, it is important to perform simulations using theoretically sound lattice quarks, such as Wilson-type quarks. Here we note that, until recently, the O(4) scaling property expected around the chiral transition of two-flavor QCD [10][11][12] had been observed only with Wilson-type quarks [13,14]. Quite recently, O(N ) scaling behavior 1 was observed with an improved staggered quark by letting the light quark mass be much lighter than the physical u and d quark masses [15]. Therefore, some of the chiral properties around the transition temperature may be easier to extract with Wilson-type quarks.
One reason for which Wilson-type quarks have not been widely used in the study of hot/density QCD is that the computational cost is larger than that for staggered-type quarks, especially at small quark masses. Previous studies with Wilson-type quarks were limited to large quark masses and the case of two-flavor QCD at vanishing chemical potentials [14,16]. The WHOT-QCD Collaboration is pushing forward studies of lattice QCD at finite T and μ using improved Wilson quarks [17] coupled to RG-improved Iwasaki glues [18,19]. We want to extend the studies to more realistic 2 + 1 flavor QCD at finite chemical potentials with physical light quarks. Towards this goal, we made a series of simulations by implementing and developing efficient methods for Wilson-type quarks. We developed the T -integration method to make the fixed-scale approach applicable [20], tested the Gaussian method to moderate the sign problem [21], and extended the histogram method by combining it with the reweighting technique to investigate the phase structure [22][23][24]. With these techniques, we have studied the thermodynamic properties of quark matter through the equation of state [25,26], heavy-quark free energies and screening masses [27][28][29], spectral functions [30,31], etc.
In Section 2, we review the techniques used in our investigation of finite-temperature QCD on the lattice. We also introduce the T -integration method to calculate the EOS in the fixed-scale approach. We then report our calculation of EOS with 2 + 1 flavors of improved Wilson quarks in Section 3. Sections 4, 5 and 6 are devoted to our study of finite-density QCD. We first discuss in Section 4 major methods to perform simulations of QCD at finite densities. We introduce our approach using the cumulant expansion in a hybrid version of the Taylor and reweighting methods. We then present in Section 5 our results for the pressure and quark number susceptibility at finite densities in twoflavor QCD, using these methods. In Section 6, we introduce another method, a histogram method, to investigate the first-order transition and its boundary. We apply the method to calculate the location of the critical point where the first-order deconfining transition of QCD in the heavy quark limit turns into a crossover as the quark masses are decreased. We then study the critical point at nonzero chemical potentials. We also present our ongoing project to study finite-density QCD with light dynamical quarks by combining the histogram method with phase-quenched simulations. Finally, in Section 7, our results for the heavy-quark free energy are summarized for zero and finite densities. A short summary is given in Section 8.

Thermodynamics of QCD on the lattice
On a lattice with a size N 3 s × N t ≡ N site , the temperature of the system is given by T = 1/(N t a), where a is the lattice spacing. To vary T , we may either vary a at a fixed N t , or vary N t at a fixed a. Let us call the former the fixed-N t approach and the latter the fixed-scale approach. 1 Presumably O(2) scaling, as expected from the symmetry of staggered-type quarks.

3/35
PTEP 2012, 01A104 S. Ejiri et al. In the fixed-N t approach, we can vary T continuously through a continuous variation of a. This is one reason that the fixed-N t approach has been widely used in many simulations.
The value of a is controlled by coupling parameters, which we denote as b. For QCD with Wilson-type quarks, b = (β, κ u , κ d , κ s , . . .). We first define the lines of constant physics (LCPs) in the coupling parameter space by the fixed dimensionless ratios of physical observables such as m π /m ρ = m π a/m ρ a. Here, to remove additional dependence on T , these observables have to be measured at T = 0. An LCP represents a physical system at different values of a. In two-flavor QCD with improved Wilson quarks, LCPs defined by m π /m ρ are determined in Refs. [16] and [27]. In 2 + 1 flavor QCD, we have to fix one more dimensionless ratio, such as m K /m ρ or m η ss /m φ . Our world corresponds to the LCP with m π /m ρ ≈ 135/770, etc.
The beta function a d b/da is defined as the variation of b along an LCP. In the fixed-N t approach, we vary T of a given physical system by varying b along an LCP on a lattice with a fixed value of N t .
The energy density and the pressure p of the system are given by derivatives of the partition function Z in terms of T and the physical volume V = (N s a) 3 : where · · · sub is the thermal average with zero temperature contribution subtracted for renormalization. To vary T and V independently, we need to introduce anisotropic lattices. When a s and a t are lattice spacings in the spatial and temporal directions, V and T are given by V = (N s a s ) 3 and T = 1/(N t a t ). Then, in principle, (2.1) can be evaluated by independent variations of a s and a t . However, this requires a systematic study of a set of physical observables on anisotropic lattices, varying both a s and a t , which is quite demanding.
Here, we note that the combination is given in terms of a uniform rescaling a t ∂ ∂a t + a s ∂ ∂a s , which can be realized without introducing anisotropic lattices. We thus obtain on isotropic lattices where S is the lattice action. The coefficient a d b/da is just the beta-function of the system, whose nonperturbative values can be determined by simulations on isotropic lattices. Equation (2.3) enables us to study − 3 p nonperturbatively without introducing anisotropic lattices. The combination − 3 p is nothing but the trace of the energy-momentum tensor, called the trace anomaly. − 3 p vanishes for free gases, but will have nontrivial values with interacting matter.

Fixed-N t approach and integration method
In order to obtain and p separately, we need one more independent input. The most widely used is the integration method [32], with which we can determine the pressure p nonperturbatively through 4 an integration in the coupling parameter space: This relation is obtained by differentiating and then integrating the thermodynamic relation p = (T /V ) ln Z in the coupling parameter space. The integration path can be freely chosen, as long as the initial point b 0 is located in the low-temperature phase such that p( b 0 ) ≈ 0. See Appendix A of Ref. [16] for a concrete demonstration of the path-independence. Several points need to be kept in mind with the fixed-N t approach, as follows: (i) When we fix N s , the spatial volume V = (N s a) 3 is varied simultaneously as we vary T . In the high T region, V may be quite small with a fixed N s . To keep V around a fixed value, N s has to be increased as we increase T . (ii) At low T , the lattice may be coarse. To ensure asymptotic scaling around the QCD transition temperature, a large N t together with improvement of the lattice action is mandatory. (iii) For the zero-temperature subtraction, we have to carry out zero-temperature simulations at all of the finite-temperature simulation points. Together with systematic zero-temperature simulations to determine the LCP over a wide range of coupling parameter space, an indispensable fraction of the total computational cost is required to carry out systematic zero-temperature simulations in the fixed-N t approach.

Fixed-scale approach and T -integration method
In the fixed-scale approach, T is varied by varying N t at a fixed b (and thus at a fixed a), i.e., the simulations are done at the same b point for all values of T . Therefore, all the simulations are automatically on an LCP without fine-tuning. Furthermore, the T = 0 subtractions can be done by a common T = 0 lattice. We may even borrow high statistic configurations at T = 0 on the International Lattice Data Grid (ILDG) [33]. Therefore, we can largely reduce the cost for the zero-temperature simulations with the fixed-scale approach [20,34].
On the other hand, the conventional integral method of obtaining p by integration in the coupling parameter space is inapplicable, because data are available at only one b point in the fixedscale approach. To overcome the problem, we have developed a new method, the T -integration method [20]. Using a thermodynamic relation at vanishing chemical potential, we obtain with p(T 0 ) ≈ 0. When we vary T by varying b along an LCP, the integral in (2.5) is equivalent to that in (2.4) with the integration path chosen to be on this LCP. However, (2.5) allows us to integrate over T without varying b.
In the fixed-scale approach, various values of T are achieved by varying N t . Because N t is discrete, we have to interpolate the data with respect to T to carry out the integration of (2.5). We need to check the magnitude of systematic errors from the interpolation of the trace anomaly in T . Application of the method to finite μ is straightforward when we reweight from μ = 0.
We find that the fixed-scale approach is complementary to the conventional fixed-N t approach in many respects, as follows. At very high temperatures, typically at T 1/3a, the fixed-scale approach suffers from lattice artifacts as discussed below, while the fixed-N t approach can keep N t finite and can keep the lattice artifact small by adopting a sufficiently large N t . On the other hand, at small T , typically at T T c , the fixed-scale approach can keep a small at larger cost due to the large N t , while the fixed-N t approach suffers from lattice artifacts due to large a.   3. Test of the fixed-scale approach armed with the T -integration method in quenched QCD [20]. Left: trace anomaly on an anisotropic lattice (a2) compared with an isotropic lattice with similar spatial lattice spacing and volume (i2). Right: energy density and pressure with the T -integral method. The shaded curves represent the results of the conventional fixed-N t method at N t = 8 [35]. See ref. [36] for a recent result.
Another attractive point of the fixed-scale approach in a study with improved Wilson quarks is that, unlike the case of the fixed-N t approach, we can keep the lattice spacing small at all temperatures and thus can avoid extrapolating the nonperturbative clover coefficient c SW to coarse lattices on which the improvement program is not quite justified.
It should be kept in mind that the fixed-scale approach is not applicable at very high temperatures where, besides the artifacts due to a quite small value of N t , the lattice spacing a may also be too coarse to resolve thermal fluctuations [20,29]. To get an idea about the latter effects, we estimate a typical length scale of thermal fluctuations by the thermal wavelength λ ∼ 1/E, where E is an average energy of massless particles at finite T . We find E ∼ 3T ζ(4)/ζ (3) ∼ 2.7T for the Bose-Einstein distribution and E ∼ 3T ζ(4)/ζ (3) × 7/6 ∼ 3.15T for the Fermi-Dirac distribution. We then obtain λ ∼ 1/3T . Thus, the data at T 1/3a for which a λ should be taken with care [29]. We have tested the fixed-scale approach and the T -integration method in quenched QCD [20]. The main results are summarized in Fig. 3. Comparing the EOS obtained on various lattices as well as the result from the fixed-N t approach on large lattices, we find that the fixed-scale approach is reliable and powerful in calculating EOS, in particular at low and intermediate temperatures. The systematic errors due to the interpolation in T are well under control in these studies. The EOS from the fixedscale approach was shown to be very consistent with that from the fixed-N t approach with large N t (N t 8), except for the high-temperature limit where the fixed-scale approach suffers from lattice discretization errors.
We adopt the fixed-scale approach to calculate the EOS in 2 + 1 flavor QCD in Section 3. We also compute heavy-quark free energies in 2 + 1 flavor QCD with the fixed-scale approach in Section 7.2.

Equation of state in 2 + 1 flavor QCD with improved Wilson quarks
A systematic study of finite-temperature QCD with improved Wilson quarks was made by the CP-PACS Collaboration around the beginning of this century for the case of two-flavor QCD at vanishing chemical potentials using the fixed-N t approach [14,16]. From a series of systematic simulations, they determined the phase structure and LCPs, confirmed the O(4) scaling, and obtained the EOS along several LCPs in the range m π /m ρ 0.65 around and above the pseudocritical temperature T pc on N t = 4 and 6 lattices. We extend the study to 2 + 1 flavor QCD. By adopting the the fixed-scale approach, we use zerotemperature configurations generated by the CP-PACS+JLQCD Collaborations [37][38][39]. Our lattice action consists of the RG-improved Iwasaki gauge action [18,19] and the clover-improved Wilson quark action [17]: 3) where c 1 = −0.331 and c 0 = 1 − 8c 1 . We set κ u = κ d ≡ κ ud and use the clover coefficient c SW nonperturbatively determined by the Schrödinger functional method in [39].
is the lattice field strength, with f x,μν the standard clover-shaped combination of gauge links. Hadronic properties with this action have been studied down to the physical point by the CP-PACS, JLQCD and PACS-CS Collaborations [37][38][39][40][41][42][43].
As the first determination of the EOS with Wilson-type quarks in 2 + 1 flavor QCD, we study at (β, κ ud , κ s ) = (2.05, 0.1356, 0.1351), which corresponds to the smallest lattice spacing and the lightest u and d quark masses (m π /m ρ 0.63) among the zero-temperature configurations generated by the CP-PACS+JLQCD Collaborations [37,38]. The s quark mass is set around its physical point (m η ss /m φ 0.74). The u and d quark masses are still much larger than their physical values. A study at the physical point [43] is reserved for the next step. The hadronic radius at this simulation point is estimated to be r 0 /a = 7.06(3) [44]. Setting the lattice scale by r 0 = 0.5 fm, we estimate 1/a 2.79 GeV (a 0.07 fm). The lattice size is 28 3 × 56 with N s a 2 fm.

Beta functions
Using the same coupling parameter values as the zero-temperature simulation, we have generated finite-temperature configurations on 32 3 × N t lattices with N t = 4, 6, . . . , 16 [25,26]. This range of N t corresponds to T 170-700 MeV.
To evaluate the trace anomaly using (2.3), we need the beta functions a(dβ/da) and a(dκ f /da) ( f = ud and s). These beta functions can be determined nonperturbatively through the coupling parameter dependence of zero-temperature observables. We use the data for am ρ , m π /m ρ and m η ss /m φ at 30 simulation points from the CP-PACS+JLQCD zero-temperature configurations [37,38] to extract the three beta functions. The first observable am ρ sets the scale. A naive method to obtain the beta functions is to fit the data of these observables as functions of the coupling parameters (β, κ ud , κ s ), and invert the matrix of the slopes ∂(am ρ )/∂β, etc. However, because the values for a(dκ f /da) are numerically much smaller than a(dβ/da), the former produce large relative errors through the matrix inversion procedure, i.e., errors for large components contaminate and dominate the errors for small components. On the other hand, from the previous experience of two-flavor QCD with improved Wilson quarks in the fixed-N t approach [16], we expect that, although the  a(dκ f /da) values are small, the quark contribution is important in the trace anomaly. Therefore, a precise determination of a(dκ f /da) is required.
To avoid the matrix inversion procedure, we instead fit the coupling parameters b = (β, κ ud , κ s ) as functions of three observables am ρ , m π /m ρ ≡ X , and m η ss /m φ ≡ Y . Consulting the overall quality of the fits, we adopt the following third-order polynomial function for the observables in this study: The global fits for each component of b are independent and have dof = 10. We find reasonable The result for β is shown in the left panel of Fig. 4. In this study, we define LCPs by m π /m ρ and m η ss /m φ at T = 0. Therefore, the beta functions are extracted as a dβ/da = (am ρ ) ∂β/∂(am ρ ), etc., in terms of the coefficients c 1 , c 2 , c 5 , c 8 , · · · in (3.5). The resulting beta functions for our LCP (m π /m ρ = 0.6337, m η ss /m φ = 0.7377) are shown in the right panel of Fig. 4 as functions of β. Their values at β = 2.05 are used to determine the trace anomaly.
As the variable to set the scale, we may alternatively use am π , am K or am K * instead of am ρ in (3.5). We use this freedom to estimate a systematic error. Taking the results from am ρ as the central value, we obtain at our simulation point, where the first brackets are for statistical errors and the second ones are for systematic errors [26]. 2 We find χ 2 /dof = 1.63, 1.08, and 1.69 for the fit of β, κ ud , and κ s , respectively. Here, χ 2 is evaluated using a standard deviation for each coupling parameter estimated by the error propagation rule using the errors of the observables and the partial derivatives of the resulting fitting function (3.5) with respect to the observables, neglecting the covariance among the observables.

Equation of state
We now calculate the trace anomaly by (2.3): where c f = 2 for f = ud and 1 for f = s. We evaluate the spatial traces in (3.8) and (3.9) by the random noise method with U(1) random numbers [21]. The number of voice vectors is 1 for each of the color and spinor indices. Our result for the trace anomaly ( − 3 p)/T 4 is shown by a thick red curve in Fig. 5. We note that the peak height of about 7 at T = 199 MeV (N t = 14) is roughly consistent with recent results for highly improved staggered quarks obtained with the fixed-N t approach at N t = 6-12 [45,46]. The shape of ( − 3 p)/T 4 suggests that T pc locates between 174 and 199 MeV.
Carrying out the T -integration (2.5) of the trace anomaly, we obtain the pressure p shown in Fig. 5. The energy density is calculated from p and − 3 p.
Apart from the larger errors, our EOS looks roughly consistent with recent results with highly improved staggered quarks near the physical point [45,46]. We note that our peak is slightly higher. This is consistent with the fact that our light quark masses are heavier than their physical values: the experience with improved staggered quarks suggests that the peak becomes slightly higher as the light quark masses are increased (see, e.g., Ref. [45]).
Our final goal is to extend the study towards the physical point, using the on-the-physical-point configurations generated by the PACS-CS Collaboration [43]. From the present study at m ud heavier than the physical point, we encountered a couple of issues, as follows. The EOS errors in Fig. 5 are larger than those obtained with the fixed-N t approach [45,46]. Besides smaller statistics, this is due to the large statistical error in ( − 3 p)/T 4 at T 200 MeV, which is caused by the enhancement factor N 4 t in (3.7) (note that S is proportional to N t ). To obtain an accurate EOS at low temperatures, we need a large statistics of O(N 7 t ) 3 . This is, however, an unavoidable step to suppress discretization errors, and the issue is common with the fixed-N t approach. Another source of systematic errors in EOS is the limited resolution in T due to the discrete variation of N t in the fixed-scale approach. Note that the lattice spacings in full QCD studies are usually coarser than those in quenched studies. Furthermore, in the present study, N t is limited to be even due to the simulation program set we have used. To improve the resolution in T , we need simulations at odd values of N t and/or a finer lattice spacing a. An alternative way would be to combine results at different a, since we can choose a to be fine with the fixed-scale approach and thus, after confirming that the discretization effects are sufficiently small in the observables under study, we may combine the results at different a to more smoothly interpolate in T . We leave these trials to a forthcoming study with lighter quarks.

Finite-density QCD on the lattice
Next, let us move on to the issues of finite-density QCD. We consider the action given by ( Here, the theory at μ f = 0 is known to have a serious problem. In a Monte Carlo simulation, we generate configurations of link variables {U x,μ } with the probability in proportion to the Boltzmann where κ = (κ u , κ d , . . .) and μ = (μ u , μ d , . . .). At μ = 0, the quark determinants are real due to the γ 5 -hermiticity of the quark matrices, However, when μ f = 0, the γ 5 -hermiticity relation changes to and thus det M f becomes complex unless μ f = 0. Because the Boltzmann weight has to be real and positive, we cannot perform a Monte Carlo simulation at μ f = 0 directly. Various methods have been proposed to study finite-density QCD avoiding the complex weight problem. However, at present, all of them are applicable in small μ f /T regions only. In the following subsections, we introduce these methods, mainly focusing on the Taylor expansion and reweighting methods we use.

Taylor expansion method
The simplest approach to studying finite-density QCD avoiding the complex weight problem is the Taylor expansion method, in which physical quantities are expanded in terms of μ f /T around μ = 0 [47][48][49][50][51]. For example, the pressure p = (T /V ) ln Z is expanded as in three-flavor QCD. The Taylor coefficients c i jk can be evaluated by a conventional Monte Carlo simulation at μ = 0 which is free from the complex weight problem. We expect that QCD in the hightemperature limit is well described by a free gas of quarks and gluons, in which p consists of terms proportional to μ 2 f and μ 4 f only. Therefore, the expansion will converge well in the high-temperature region.
Other observables can also be calculated similarly. The quark number density n f is given by where T is fixed in the differentiations. When we define the light quark number density as n q = n u + n d , the susceptibilities of the light quark number density (χ q ) and the strange quark number density (χ s ) are given by respectively. The susceptibility of the isospin number can also be given as These quantities are expanded around μ = 0 in terms of c i jk defined in (4.5). The trace anomaly is given by Eq. (2.3). The entropy density is given by The chiral condensate is defined by the derivative of ln Z with respect to the quark mass. Their Taylor expansion can also be derived.

Reweighting method and sign problem
Another popular approach to finite-density QCD is the reweighting method [52][53][54][55][56][57], using the reweighting technique [58][59][60] in finite-density QCD. Using the configurations generated at μ = 0, The denominator is the ratio of the partition functions at finite μ and μ = 0, has a complex phase e iθ . When the fluctuation ofθ is larger than O(π/2), a reliable calculation of the numerator and denominator in Equation (4.9) is difficult. This difficulty is called the "sign problem (complex phase problem)". We actually encounter large fluctuations ofθ at large μ f and/or large lattice volume. It is worth rewriting the denominator of Equation (4.9) in terms of the distribution function for P ≡ −S g /(6N site β) (which is the plaquette if S g is the standard gauge action) and the logarithm of the absolute value of the reweighting factorF ≡ ln where w 0 is the distribution function for the phase-quenched system, and e iθ P,F is the expectation value of the operator e iθ at μ = 0 with fixed values ofP andF to P and F: (4.13) By measuring the histogram ofP andF in a phase-quenched simulation, we can determine w 0 around the peak of the histogram. However, in the calculation of Equation (4.11), precise information on w 0 is required around the maximum of the integrant, which may deviate from the peak of w 0 due to the factor e iθ P,F e F . To calculate w 0 in a wider range of (P, F), we combine the results of phasequenched simulations at different points in the coupling parameter space, adopting the reweighting formulae for the phase-quenched theory [21]. Further demonstration of such calculations will be given in Section 6.
For Equation (4.11), we also need to estimate e iθ P,F . Because the total distribution function is real and positive in finite-density QCD due to the charge conjugation symmetry, the imaginary part of e iθ P,F must be averaged out when the statistics is infinite. Since the imaginary part is the source 12  of the sign problem, we may remedy or mitigate the problem if we could find a method in which the imaginary part is removed and the real part is reliably estimated. In the next subsection, we show that it is useful to consider a cumulant expansion of e iθ P,F in which ln e iθ P,F is separated into real and imaginary parts [21,61].

Cumulant expansion and Gaussian approximation
For simplicity, let us consider the case of N f degenerate quarks in the following. A crucial step in handling the fluctuations in the phaseθ is to introduce an appropriate definition ofθ removing the ambiguity of a complex phase with modulo 2π . We uniquely define the complex phase by the Taylor expansion asθ where the derivatives of ln det M can be unambiguously expressed in terms of M −1 and derivatives of M. Note that the expectation value ofθ defined by Equation (4.14) is not restricted to the range (−π, π), and the maximum value of |θ | is proportional to the volume of the system. The conventional complex phase in the range (−π, π) is recovered by taking the principal value ofθ with modulo 2π . Typical results for the histogram ofθ are shown in Fig. 6 for N f = 2 QCD with improved Wilson quarks [21], where the power series in (4.14) is evaluated up to (μ/T ) 3 , discarding O (μ/T ) 5 terms 4 . We see that the width of the distribution becomes wider as μ/T is increased, indicating a more severe sign problem.
To mitigate the sign problem, we evaluate e iθ P,F by the cumulant expansion [61]: where θ n c is the n th -order cumulant: A key observation in handling the cumulant expansion is that θ n c = 0 for any odd n due to the symmetry underθ → −θ. This implies that, provided that the cumulant expansion converges, e iθ P,F is guaranteed to be real and positive and the sign problem is resolved. The sign problem is transformed into the convergence problem of the cumulant expansion in this approach.
As shown in Fig. 6, we find that the distribution ofθ is well described by a Gaussian function up to μ/T ∼ O (1). The Gaussian distribution ofθ is also observed with improved staggered quarks in two-flavor QCD [61], and was discussed in Refs. [62,63] too. The Gaussian distribution means that the cumulant expansion (4.15) is dominated by the lowest nontrivial order of θ 2 , and thus the expansion is well converged. Corrections to the Gaussian distribution can be incorporated by taking higher-order terms in the cumulant expansion.
Here, we note thatθ = O(μ/T ) from Equation (4.14). Therefore, if we take into account the cumulants up to the n th order, the truncation error does not affect the Taylor expansion up to O ((μ/T ) n ). This means that the convergence of the cumulant expansion is closely related to the convergence of the Taylor expansion in μ. The Gaussian approximation is valid, at least at small μ, and higher-order cumulants may become visible at large μ. The applicability range of the Gaussian approximation in μ has to be checked for each case by calculating the higher-order cumulants.
We now argue that the range of applicability does not change with the system volume on sufficiently large lattices when the correlation length of the system is finite [21]. The expansion coefficients forθ in Equation (4.14) are given by traces of products of M −1 and ∂ n M/∂(μ/T ) n over the spatial positions. For example, the first coefficient is given by the trace of N f [M −1 (∂ M/∂(μ/T ))]. We note that the diagonal elements of this matrix are the local quark-number density operators [∼ψ(x)γ 0 ψ(x)] at μ = 0. When the correlation length of the local operators is finite and is much shorter than the system size, we may decompose the trace into a summation of independent contributions from spatially separated regions. The same discussion is applicable to higher-order coefficients too. Then, the phaseθ may be written as a sum of local contributions from spatially separated regions, which are statistically independent of each other, i.e.θ ≈ xθ x , whereθ x is the contribution from a spatial region labeled by x. The average of exp(iθ) is thus This suggests that all cumulants θ n c ≈ x θ n x c increase in proportion to the volume as the volume increases, that is, θ n c ∝ volume for any n, in contrast with a naïve expectation thatθ n may be proportional to (volume) n sinceθ is proportional to the volume. 5 Therefore, while the width of the distribution, i.e. the phase fluctuation, increases in proportion to the volume, the ratios of the cumulants are independent of the volume-the higher-order terms in the cumulant expansion are well under control in the large volume limit. 5 This property of the higher-order cumulants can also be understood from the effective potential V eff (P, F) = − ln w 0 (P, F) − ln e iθ (P,F) , which will be studied in Section 6. Because V eff and ln w 0 are proportional to the system volume, each term in the expansion of ln e iθ (P,F) = n (i n /n!) θ n c should not increase faster than (volume) 1 for any n. Otherwise, V eff becomes singular in the thermodynamic limit. Thus, the range where the cumulant expansion is applicable is independent of the volume. This is good news in attacking the sign problem, which is known to become more severe with increasing lattice volume. Furthermore, we find that, when the system size is sufficiently larger than the correlation length, the distribution ofθ/volume tends to a Gaussian distribution, since the distribution function of an average over many independent variablesθ x tends to a Gaussian function by the central limit theorem.

Other approaches
Besides the Taylor expansion method and the reweighting method, as well as combinations thereof, various methods have been proposed as alternative approaches to study QCD at finite densities. One approach is to perform analytic continuation from simulations at imaginary chemical potentials [64][65][66][67][68]. For a complex μ f , Equation Therefore, when μ f is purely imaginary, the Boltzmann weight is real, and we can simulate the system without the sign problem. Using results of simulations performed at imaginary μ, information at small real μ can be obtained by an analytic continuation. The analytic continuation is usually based on a Taylor expansion in terms of μ around μ = 0, i.e. we fit observables obtained at imaginary μ with the Taylor expansion ansatz and extrapolate the resulting fitting function to a small real μ. Improvements of the analytic continuation procedure have also been discussed in Refs. [69][70][71][72] to obtain results in a wider range of real μ. Systematic high precision simulations over a wide imaginary-μ range are required for a reliable analytic continuation.
Another approach is to construct the canonical partition function Z C (T, N ) by fixing the total quark number N or the quark number density [55][56][57][73][74][75][76][77][78]. Using the canonical partition function, we can also discuss the effective potential as a function of the quark number. The relation between the grand canonical partition function Z (T, μ) and the canonical partition function Z C (T, N ) is given by for the degenerate N f -flavor case. Because this is a Laplace transformation, Z C is obtained from Z by an inverse Laplace transformation. To find N that makes the largest contribution to Z , it is worth introducing an effective potential as a function of N , where f is the Helmholtz free energy. V eff is useful to study the nature of phase transitions [78].

Two flavors of improved Wilson quarks at finite density
Using the Taylor expansion method and the Gaussian method discussed in the previous section, we made the first study of finite-density QCD with Wilson-type quarks [21]. We study two-flavor QCD with RG-improved gauge action and clover-improved Wilson quark action. A systematic study of finite-temperature QCD with this action has been made by the CP-PACS Collaboration [14,16]. The phase diagram at μ = 0 as well as LCPs determined with m π /m ρ are given in Refs. [14,16,21]. We extend the study to μ = 0.

Taylor expansion for EOS
We study the pressure and quark number densities defined by Eqs. potential μ I = (μ u − μ d )/2, the quark number and isospin susceptibilities are given by which measure the fluctuations in baryon and isospin numbers in the medium [79].
We study the isosymmetric case μ u = μ d = μ, μ I = 0. The pressure is given by Here, c 0 (T ) is the pressure at μ = 0 and has been computed by the CP-PACS Collaboration [14,16].
The susceptibilities are then expanded as where In this study, we compute the Taylor expansion coefficients for the second and fourth derivatives in (5.2). This enables us to compute χ q and χ I to the lowest nontrivial order in μ.

Random noise method
To evaluate the Taylor coefficients (5.2) and (5.4), we calculate D n = N f ∂ n ln det M ∂(μa) n μ=0 (5.5) up to n = 4 [21]. We thus study In terms of D n , we have . . . Note that D n is real for even n and purely imaginary for odd n [47,48].
These traces, say trA, can be evaluated by the "random noise method". In this method, instead of calculating the diagonal elements A ii individually for all i, we calculate η † Aη for several random noise vectors η. The contributions of the off-diagonal elements A i j (i = j) in these quantities are removed by averaging over the random noises, η * i η j = δ i j . As is clear from the construction, the random noise method is effective when the contaminations from off-diagonal elements are small. Because the propagator (M −1 ) x,y decreases rapidly with increasing |x − y|, the random noise method will work well to suppress small contaminations of spatially off-diagonal elements. On the other hand, the off-diagonal elements in the color and spinor 16 indices are from the same spatial point and thus are not suppressed by |x − y|-they are suppressed only by 1/ √ N noise , where N noise is the number of noise vectors. This motivates us to apply the random noise method for the spatial coordinates only. For the trace over the color and spinor indices, we just repeat the calculation generating the random noise vectors for each color-spin index 6 .
For a product of traces, the random noise vectors for each trace must be independent. We compute such a product by subtracting the contribution of the same noise vector from the naive product of two noise averages for each trace. This effectively increases the number of noises to N noise · (N noise − 1) for the products and thus suppresses their errors due to the noise method.
In practice, N noise can be small when the error due to the smallness of N noise is smaller than the statistical errors from the averaging over the configurations. The required number of noise vectors depends on each operator. In Ref. [21], we have tested the random noise method in the evaluation of D n (with n = 1-4). We find that D 1 has larger fluctuations under noise vector variation than D 2 etc., and the error in

Quark number densities and susceptibilities at μ = 0
Simulations are done on a 16 3 × 4 lattice in the fixed-N t approach in the range β = 1.50-2.40 (T /T pc ≈ 0.8-3 or 4 on LCPs corresponding to m π /m ρ = 0.65 and 0.80, where T pc is the pseudocritical temperature for each LCP. See Ref. [21] for details. Our results with the standard Taylor expansion method are shown in Fig. 7 for the LCP at m π /m ρ = 0.65. The left panel is for the finite density correction p/T 4 ≡ p(μ)/T 4 − p(0)/T 4 of the pressure at finite μ. T 0 is T pc at μ = 0. Recall that we have a crossover at T pc at μ = 0. The pressure changes more sharply as μ is increased. When we increase μ, p/T 4 becomes the same size as p/T 4 at μ = 0 around μ/T ∼ O(1). In the right panel of Fig. 7, we show the results of the quark number susceptibility χ q (μ) at μ = 0. We see that a peak seems to be formed when we increase μ. However, in spite of various improvements in random noise estimators etc., as discussed in the previous subsection, the statistical errors due to the complex phase fluctuation of the quark determinant are still a bit too large to draw a definite conclusion about the peak.
In order to suppress the errors from the phase fluctuation, we apply the cumulant expansion method discussed in Section 4.3. We compute the quark determinant ratio in Eq. (4.11) by the Taylor expansion up to O(μ 4 ), 6 Because a staggered-type quark does not have a spinor index at a spatial point, the number of off-diagonal elements is only 6 in the color 3 × 3 matrix, and the contamination of off-diagonal elements is less serious. This is one reason that the random noise method is used more naively with staggered-type quarks. However, with Wilson-type quarks, the number of off-diagonal elements with similar magnitude in the quark matrix is 11 times larger than the diagonal one. Therefore, the color-spinor index should be treated more carefully with Wilson-type quarks.   with N max = 2, and estimate the phase factor by the second-order cumulant assuming the Gaussian distribution of θ , e iθ ≈ exp[− θ 2 /2]. We moreover shift β from the simulation point β 0 such that the statistical error due to the F-integration is minimized in Eq. (4.11). We perform a fit of the resulting pressure (= free energy) in terms of μ to obtain χ q . The results for p/T 4 and χ q /T 2 are shown in Fig. 8. We find that the statistical errors are appreciably suppressed by these methods. We also note that, although the simulations at different T are independent, the T -dependences of p and χ q are smooth and natural. We can now clearly identify a sharp peak in χ q /T 2 that appears around T pc when μ/T ≥ O(1). The peak becomes higher as μ increases. These points are consistent with observations with staggered-type quarks and suggest a critical point at finite μ.
In contrast with the peak in χ q , the isospin susceptibility shows no sharp peak, as shown in the left panel of Fig. 9. This is in accordance with the expectation that χ I is analytic at the critical point since the iso-triplet mesons remain massive.
The results at m π /m ρ = 0.80 are similar, but with lower peaks in χ q than those in Fig. 8, as shown in the right panel of Fig. 9. This may be explained in part by the expectation that the critical point locates at larger μ because the quark mass is larger than that for m π /m ρ = 0.65. See Ref. [21] for more discussions.   9. Results of two-flavor QCD at finite density [21]. Left: Isospin susceptibility χ I at m π /m ρ = 0.65 by the standard Taylor expansion method. Right: Quark number susceptibility χ q at m π /m ρ = 0.80 by the combined hybrid and Gaussian methods.

Histogram method and QCD phase structure at zero and finite densities
In studying the QCD phase structure shown in Figs. 1 and 2, identification of first-order transition region is quite important. Among several methods to study the nature of phase transitions, the probability distribution of physical observables provides us with the most intuitive. The probability distribution function is defined as the generation rate of configurations with fixed expectation values of physical observables. Double or multiple peaks in the probability distribution function for observables that are sensitive to the phase, such as the energy density, chiral order parameter, etc., give a signal of first-order phase transition. One problem is that, in order to trace the variation of the shape of a probability distribution function, we need statistically reliable data on the distribution function over a wide range of expectation values. In the identification of a first-order transition, correct evaluation of the double-peak distribution requires quite a long simulation with sufficiently many flip-flops among different phases [see, e.g., Ref. [80]]. This is computationally quite demanding with dynamical quarks.
Here, we note that calculation of the probability distribution function is also required for the reweighting method for the components of the action [58][59][60]. Therefore, when we adopt observables that are the components of the action, the computation of the distribution function at different points in the coupling parameter space is straightforward with the reweighting method [61]. This helps us to obtain the probability distribution function over a wide range of expectation values. Therefore, the method is quite powerful for studying the location of first-order transitions.
We apply the method to explore the phase structure of QCD. In this section, we introduce the method, which may be viewed as a variant of the histogram method or the density of state method [81][82][83][84], and test it in the heavy-quark region of QCD [22,23]. We then present our ongoing project to study finite-density QCD with light dynamical quarks by combining the histogram method with phase-quenched simulations [24].

Histogram method
As a demonstration, let us consider the simplest lattice QCD: the combination of plaquette gauge action with unimproved Wilson quarks. 19 2) where N site = N 3 s × N t is the lattice volume and is the plaquette. Note that M does not depend on β. 7 Denoting the values of the operators (P, . . .) as (P, . . .), the probability distribution function for (P, . . .) is defined by (6.4) where "· · · " in the r.h.s. means the product of delta functions for other operators. We now define the effective potential as V eff (P, . . . ; β, κ, μ) = − ln w(P, . . . ; β, κ, μ). (6.5) Note that P represents the freedom of gauge internal energy, and thus should be sensitive to the phase structure of the system. A useful property of the plaquette distribution function and the effective potential is w(P, . . . ; β, κ, μ) = w(P, . . . ; β 0 , κ, μ) e 6(β−β 0 )N site P , (6.6) We thus find that The κ and μ-dependences of V eff can also be computed by the reweighting method as follows: R(P, . . . ; κ, μ, κ 0 , μ 0 ), (6.9) where the reweighting factor R is evaluated as Note that R is independent of β and thus can be evaluated at any β. By adjusting and combining β, we can obtain precise values of R in a wide range of P, . . . .

QCD in the heavy-quark region
We first test the method in the heavy-quark region. As shown in Fig. 2, we have the first-order deconfinement transition of the SU(3) Yang-Mills theory in the limit of infinite quark masses. The transition is expected to turn into a crossover when we decrease the quark masses. We study the boundary of the first-order region by the histogram method. To take advantage of light computational costs in quenched simulations, we choose κ 0 f = μ 0 f = 0 ( f = 1, . . . , N f ). For simplicity, let us consider the case of degenerate quarks: κ f = κ, μ f = μ ( f = 1, . . . , N f ). Generalization to nondegenerate cases is easy.
To the lowest order of the hopping parameter expansion, the quark determinant ratio appearing in the r.h.s. of (6.10) is evaluated as Tr U x,4 U x+4,4 · · · U x+(N t −1)4,4 (6.12) is the Polyakov loop, andˆ R = Reˆ andˆ I = Imˆ are its real and imaginary parts. We note that the contribution ofP in the r.h.s. can be absorbed by a shift of the gauge coupling β → β + 48N f κ 4 . The term proportional toˆ I induces the complex phase at μ = 0, which is the origin of the sign problem at large μ.

Results at μ = 0 in the heavy-quark region
We are now ready to calculate the effective potential. Let us first study the case μ = 0 [22,23]. In this case, the complex phase term is absent in (6.11). 21 Fig. 10. Effective potential at μ = 0 as a function of P in the heavy-quark region [22]. Left: dV eff (P; β, 0)/d P in the heavy quark limit. Data obtained at five different values of β in the range 5.68-5.70 are converted to β = 5.69 by (6.8). Right: V eff (P; β, κ) for κ > 0, where β is adjusted such that the two minima of V eff have the same depth and the constant term of V eff is adjusted such that V eff = 0 at the central peak point P peak .
The double-well nature of the effective potential is clearly seen when we consider the effective potential for P only: Our results for dV eff (P; β, 0)/d P at κ = 0 are shown in the left panel of Fig. 10. Using (6.8), we shift the results obtained at five different β values to β = 5.69. With different β, the range of P in which we can reliably obtain V eff is different. We find that the results of dV eff (P; β, 0)/d P at different β are smoothly connected with each other by (6.8). We can thus obtain accurate values of dV eff (P; β, 0)/d P in a wide range of P.
Similarly, we calculate dV eff (P; β, κ)/d P at κ > 0 by using the reweighting factor R and (6.11). We then integrate dV eff (P; β, κ)/d P in P to get V eff (P; β, κ) as shown in the right panel of Fig. 10. We find that the double-well structure of V eff (P) becomes weaker with increasing κ, and eventually disappears at finite κ, say κ cp . Examining the shape of V eff more closely, we obtain κ cp = 0.0658(3)( +4 −11 ) for two-flavor QCD in the lowest order of the hopping parameter expansion on an N t = 4 lattice [22].
The argument can be easily generalized to the case of nondegenerate quark masses. Our results for the critical point κ cp in 2 + 1 flavor QCD are shown in Fig. 11 [22]. The top-right corner of Fig. 2 corresponds to this plot rotated by 180 • . Our results are consistent with those obtained in an effective model [86] and with a recent study using the hopping parameter expansion [87].
The expression (6.11) leads us to adoptˆ R as an additional operator for the effective potential V eff . Then,ˆ R in the r.h.s. of (6.11) is simply replaced by its expectation value R , and the reweighting factor R is just a given function of P and R in this case. We have where V 0 is the effective potential in the heavy quark limit (SU(3) pure gauge theory). The argument β in ∂ V eff /∂ R and ∂ V 0 /∂ R is omitted in (6.14) since they are independent of β due to (6.7). Note that, besides known overall constants [the last terms in (6.13) and (6.14)], the dependences of ∂ V eff /∂ P and ∂ V eff /∂ R on P and R are independent of β and κ.
Because R represents the freedom of heavy-quark free energy, we expect it be sensitive to the phase structure of the system. Around the first-order transition point, we will have a double-well structure of V eff in the 2D plane (P, R ). To study the phase structure, it is useful to examine the curves ∂ V eff /∂ P = 0 and ∂ V eff /∂ R = 0. From (6.13) and (6.14), these curves at different (β, κ) correspond to different contour curves of ∂ V 0 /∂ P and ∂ V 0 /∂ R . A contour plot for ∂ V 0 /∂ P and ∂ V 0 /∂ R is given in Fig. 12. When the curves ∂ V eff /∂ P = 0 and ∂ V eff /∂ R = 0 cross at only one point, we have just one minimum of V eff . In this case, we have no first-order transitions around this (β, κ). On the other hand, when we have three intersection points, we have two minima and one saddle point, implying the existence of the first-order transition. In particular, from the merger of three intersection points, we can determine the critical point where the first-order transition line terminates. From Fig. 12, we find that the S-shaped curve of ∂ V eff /∂ R = 0 leads to three intersection points 23 at small κ, and that the S-shape becomes weaker with increasing κ. A preliminary estimate for the critical point is κ cp ≈ 0.0690(7) [23], shown by the thick contour curves in Fig. 12. We are currently testing a refinement of the method to extract smoother contour curves.
To examine the quality of the hopping parameter expansion, we have studied the effect of the nextleading order terms in the evaluation of the critical point κ cp (WHOT-QCD Collaboration, manuscript in preparation and Ref. [85]). To this order, we need to incorporate κ 6 -loops with length six and generalized Polyakov loops with length N t + 2 in the quark determinant ratio. The effects of κ 6 -loops may be absorbed by a shift of β. Examining the effects of generalized Polyakov loops, we find that κ cp only shifts by about 3% on the N t = 4 lattice due to the next-leading order terms. We also find that the contribution of the next-leading order terms becomes comparable to that of the leading order terms at κ ∼ 0.18. Because κ cp is much smaller than this, we conclude that the hopping parameter expansion is well valid up to κ ∼ κ cp . Accordingly, an estimation of the pseudoscalar meson mass around κ cp leads to m π ≈ T /0.023 for N f = 2 [22], i.e., m π ∼ 7-9 GeV with T ∼ T pc ∼ 160-200 MeV. Thus, κ cp locates well in the heavy-quark region.

Results at μ = 0 in the heavy-quark region
In order to extend the study to finite densities, we have to calculate the reweighting factor due to the complex phase, where · · · P, R is the expectation value at fixed P and R in quenched QCD. Whenθ fluctuates a lot at large μ, a reliable estimation of e iθ P, R becomes difficult (the sign problem). 8 Before evaluating e iθ P, R , let us consider the case of phase-quenched finite-density QCD, in which the complex phase term is removed in the quark determinant. In two-flavor QCD, this corresponds to the case of the isospin chemical potential, μ u = −μ d ≡ μ I . From (6.11), we find that, after shifting β → β + 48N f κ 4 , the effects of μ I are just to further shift κ → κ cosh 1/N t (μ I /T ). Therefore, we have κ I cp (μ I ) = κ cp (0)/ cosh 1/N t (μ I /T ) (6.16) for the critical point in the phase-quenched QCD to the lowest order of the hopping parameter expansion, where κ cp (0) is the critical point at μ u = μ d = 0. Note that, with increasing μ I , the critical point approaches κ = 0, where the hopping parameter expansion is reliable. We now compute the effect of the phase. To evaluate e iθ P, R , we adopt the Gaussian approximation with the cumulant expansion [21,61] discussed in Section 4.3. In the heavy-quark region, we study θ 2n We find that θ 2 c θ 4 c around the critical point [23]. This confirms the validity of the Gaussian approximation. We thus have S. Ejiri et al. When the last term in (6.18) modifies the S-shape of the curve ∂ V eff /∂ R = 0 shown in Fig. 12, κ cp (μ) deviates from κ I cp (μ). By evaluating these quantities, however, we find that the contribution of the last term in (6.18) is at most about 3% of the second term around κ I cp (μ), even in the large μ limit. Therefore, κ cp (μ) is indistinguishable from κ I cp (μ) within the current statistical errors [23]: Generalization of this result to nondegenerate cases such as the 2 Critical surfaces for the symmetric case μ u = μ d = μ s ≡ μ and a more realistic case of μ u = μ d ≡ μ ud and μ s = 0, which may be realized in heavy ion collisions, are shown in Fig. 13.

Phase-quenched simulations towards the physical point
Our challenge is to apply the histogram method to explore the phase structure of finite-density QCD in the light quark region. For the sake of notational simplicity, we consider the degenerate case in this section too. Generalization to nondegenerate cases is straightforward. When quarks are light, the Polyakov loopˆ no longer plays a decisive role in the dynamics of the system. We thus consider det M itself as the additional operator for the effective potential. det M represents the freedom of quark internal energy, and thus should be sensitive to the phase structure of the system. We denote the absolute value and the complex phase of the quark determinant as follows: where w 0 and V 0 = − ln w 0 are the distribution function and the effective potential for the phasequenched system, and the expectation value e iθ (0:β,κ,μ);P,F is evaluated with fixed P and F in the phase-quenched simulation.
In the following, let us consider a simpler case where the quark matrix M can be treated as independent of β. For example, in a study of the phase structure at a particular value of β, we may treat β in M as being fixed at that value. Physical properties such as the phase structure will not be affected by this procedure 9 . Then, e iθ (0:β,κ,μ);P,F becomes independent of β and can be evaluated at any β to cover a wider range of P and F.
Because the phase-quenched simulations in two-flavor QCD correspond to the case of isospin chemical potentials, a comment is in order about the influence of the pion-condensed phase that exists at large isospin chemical potentials [88]. In the pion-condensed phase, e iθ is expected to vanish by model studies [89,90]. According to (6.23), this means that the configurations in the pioncondensed phase make no contributions to the physics of phase-unquenched QCD-w and V eff are dominated by phase-quenched configurations out of the pion-condensed phase, and we only need to generate configurations outside the pion-condensed phase. We test the method in two-flavor QCD using the RG-improved Iwasaki gauge action (3.2) and the clover-improved Wilson quark action (3.3). According to the footnote in Section 6.1,P is identified aŝ x,μν (6.25) in the present case. We perform simulations at m π /m ρ ≈ 0.8 on an 8 3 × 4 lattice with a nonperturbatively estimated c SW . In the present study of the phase structure at each β, we treat c SW as a constant independent of β.
In the left panel of Fig. 14, we show our results of V 0 as a function of F (after integrating out P). The results of simulations at three values of μ 0 /T are translated to μ/T = 2.0 using a reweighting formula for V 0 in the μ-direction: When we calculate observables as functions of β, we need to incorporate the effects from the β-dependence in M. See [91] as a trial to incorporate the β dependence of c SW in the reweighting factor.
where R 0 is not dependent on β. From this figure, we confirm that the data obtained at different simulation points form a smooth V 0 . We can thus obtain precise V 0 over a wide range of F.
To calculate e iθ (0:κ,μ);P,F , we adopt the cumulant expansion method. A typical result for the distribution of θ is shown in the right panel of Fig. 14. We find that the distribution is well approximated by a Gaussian function. We also note that the second-order cumulant increases with increasing μ, while the fourth-order cumulant is consistent with zero within the statistical error, though the statistical error increases with μ [24]. Therefore, the cumulant expansion is well controlled by the leading term and we may reliably evaluate the complex phase factor even at these relatively high values of μ. A project towards clarification of the phase structure at the physical point is currently underway.

Heavy-quark free energy
Finally, we study the heavy-quark free energies and screening masses in QGP. The free energies for a static quark-antiquark and two static quarks characterize the inter-quark interactions in QGP, and their Debye screening masses describe the thermal fluctuation of quarks and gluons in QGP. In a phenomenological model, they are relevant to the fate of heavy-quark bound states such as J/ψ and ϒ in QGP created in relativistic heavy-ion collisions at RHIC and LHC [92,93]. On the lattice, studies in quenched QCD [94][95][96][97][98][99] and in full QCD with staggered-type quark actions [100][101][102] or Wilson-type quark actions [21,[27][28][29]103] have been made. Comparisons with analytic studies [104,105] have also been attempted.
Heavy-quark free energies on the lattice are defined through the correlation functions of the local 4 . Note that the trace over the color indices is not taken inˆ (x). With gauge-fixing, we define the free energy F R in various color channels R [106,107]:  Trˆ (x)ˆ (y) , (7.4) where r = |x − y|, and R = 1 for the color singlet Q Q channel, 8 for the color octet Q Q channel, 3 * for the color antitriplet Q Q channel, and 6 for the color sextet Q Q channel, respectively. To preserve the free energy interpretation of F R by transfer matrix theory, the gauge-fixing procedure should not include temporal links. We thus adopt the Coulomb gauge.
Above T pc , we also introduce normalized free energies V R , whose constant terms are adjusted such that they vanish at large distances r → ∞. This is equivalent to defining the free energies by dividing the r.h.s. of Equations. (7.1)-(7.4) by Trˆ (x) 2 .

Heavy-quark free energies in two-flavor QCD
We first compute the free energies (7.1)-(7.4) in two-flavor QCD at zero and finite densities [27]. We consider the case μ u = μ d = μ. We use the gauge configurations generated for the studies discussed in Section 5. As mentioned in Section 5.3, the configurations were generated on a 16 3 × 4 lattice on LCPs corresponding to m π /m ρ = 0.65 and 0.80. We thus adopt the fixed-N t approach for this study.
7.1.1. Heavy-quark free energies at μ = 0 in two-flavor QCD. The heavy-quark free energies at μ = 0 are shown in Fig. 15 for m π /m ρ = 0.65 and T ≥ T pc . The results for m π /m ρ = 0.80 are similar. We find that the inter-quark interaction is "attractive" in the color singlet and antitriplet channels and is "repulsive" in the color octet and sextet channels. We also see that, irrespective of the channels, the inter-quark interaction rapidly becomes weak at long distances as T increases, as expected from the Debye screening at high temperatures.
We find that the heavy-quark free energies in the high-temperature phase are well described by the screened Coulomb form, where α eff (T ) and m D (T ) are the effective running coupling and Debye screening mass, respectively. From the fits, we note that the color-channel dependence of the free energies can be absorbed by the 28   kinematical Casimir factor inspired by the lowest-order perturbation theory: This Casimir dominance at high temperatures has also been reported in quenched studies using the Lorentz gauge [95,96]. With the Casimir factors, α eff (T ) and m D (T ) are universal to all color channels, as shown in Fig. 16. The magnitude and the T -dependence of the Debye mass are also consistent with the next-to-leading-order thermal perturbation theory [27].

7.1.2.
Heavy-quark free energies at μ = 0 in two-flavor QCD. The Taylor expansion of V R with respect to μ q /T is given by where concrete forms of the Taylor coefficients v R n are given in Ref. [21]. The color singlet and octet channels do not have odd orders in the Taylor expansion since the QQ free energies are invariant under the charge conjugation. v R 0 is the normalized free energy at μ = 0 shown in Fig. 15. The results for the lowest nontrivial order are shown in Fig. 17. See Ref. [21] for higher orders as well as those at m π /m ρ = 0.80. From these figures, we find that, both around T pc and at higher temperatures, the sign of v R 1 is the same as that of v R 0 , whereas the sign of v R 2 is the 29/35 PTEP 2012, 01A104 S. Ejiri et al. Fig. 18. Free energies of static quarks in 2 + 1 flavor QCD at finite temperatures with the fixed-scale approach [29]. The scale was set by the Sommer scale r 0 = 0.5 fm. Left: Bare free energy in the color singlet channel. The static quark-antiquark potential V (r ) at T = 0 has been calculated by the CP-PACS and JLQCD Collaborations [38]. The fit result of V (r ) by the Coulomb + linear form is shown by the dashed gray curve. The arrows on the right side denote twice the single-quark free energy 2F Q . Right: Normalized free energies for color singlet and octet QQ channels.
Because v R 1 is absent for QQ free energies, this means that, in the leading-order of μ q , the interquark interaction between Q andQ becomes weak at finite μ q , while that between Q and Q becomes strong. In other words, QQ (Q Q) free energies are screened (antiscreened) by the internal quarks induced at finite μ q .
Taylor expansion coefficients for α eff (T, μ) and m D (T, μ) can be computed similarly. We find that the leading correction in m D (T, μ) due to finite μ is larger than the prediction from leading-order thermal perturbation theory [21].

Heavy-quark free energies in 2 + 1 flavor QCD
We now extend the study to the more realistic 2 + 1 flavor QCD [29,44]. As discussed in Section 3, we adopt the fixed-scale approach for 2 + 1 flavor QCD. We thus vary T by varying N t with fixed coupling parameters. We use the finite-temperature configurations generated in Section 3 to compute the heavy-quark free energies at m π /m ρ 0.63 and m η ss /m φ 0.74. In this study of 2 + 1 flavor QCD, we restrict ourselves to the case μ = 0.
A good feature of the fixed-scale approach is that the renormalization factors, which are determined on a zero-temperature lattice depending on the coupling parameters, are common to all T in the fixed-scale approach, because the coupling parameters are fixed for all T . This is also the case for the heavy-quark free energies. Therefore, we can directly compare the bare free energies at different T in the fixed-scale approach.
Our results for the bare free energies are shown in the left panel of Fig. 18. Plotted are the data for the color singlet QQ channel at various T in the high-temperature phase, together with the zerotemperature static quark-antiquark potential V (r ) defined by the Wilson loop operator: (7.9) We find that singlet free energy F 1 (r, T ) at any T converges to V (r ) at short distances. This is in accordance with the expectation that the short-distance physics is insensitive to temperature. With the fixed-scale approach, we directly confirm that the theoretical expectation is actually satisfied on the lattice. 10 At large distances, F 1 (r, T ) departs from V (r ) and eventually becomes flat due to Debye screening. On the right edge of the left panel of Fig. 18, we show twice the single-quark free energy defined by 2F Q = −2T ln Tr (x) at each T by arrows of the same color. We confirm that F 1 (r, T ) converges to 2F Q (T ) quite accurately.
By subtracting 2F Q , we obtain the normalized free energies shown in the right panel of Fig. 18 for the QQ channels. See [29] for the results for the Q Q channels. Performing fits with the screened Coulomb form (7.5), we confirm the Casimir dominance (7.6), as in the case of two-flavor QCD.

Gauge-independent screening masses
The Debye screening masses and the effective couplings computed in the previous subsections are dependent on the choice of the gauge. Therefore, their physical meanings are not quite clear. In Ref. [28], we have proposed a gauge-independent definition of screening masses for electric and magnetic channels.
Under the Euclidean time-reflection R and the charge conjugation C, the gluon fields transform as We call an operator magnetic (electric) if it is even (odd) under R. It is natural to extract magnetic and electric properties by decomposing observables using these symmetries [110]. Under these transformations, the Polyakov line operatorˆ (x) transforms aŝ A magnetic (electric) Polyakov line operator can be defined as the R-even (R-odd) part ofˆ (x), which can be further decomposed into C-even and odd parts aŝ 10 In the case of the conventional fixed-N t approach, a different renormalization factor is required at each T . In early studies, insensitivity at short distances was just assumed and used to adjust the constant term of F 1 (r, T ) at different T [97,100]. In more recent studies, the renormalization factors are computed by a series of zero-temperature simulations at the same coupling parameters as the finite temperature simulations. See, e.g., Ref. [108,109].   19. Comparison [28] of the screening ratio, m E− /m M+ , with predictions in the dimensionally-reduced effective field theory (3D-EFT) [111] and N = 4 supersymmetric Yang-Mills theory [112].
where r = |x − y| and Trˆ is real due to the C symmetry. Note that the conventional gaugeinvariant Polyakov loop correlation function is given by Using the configurations generated for [27], we computed these screening masss in two-flavor QCD in the high-temperature phase. We find that (i) C M+ and C E− have opposite signs, and (ii) C M+ has larger magnitude and longer range than C E− at long distances. The latter implies that the the conventional C is dominated by the magnetic sector at long distances, and thus m (T ) m M+ (T ). We also find that m M+ (T ) < m E− (T ). A comparison with dimensionally-reduced effective field theory [111] and N = 4 supersymmetric Yang-Mills theory with anti de Sitter/conformal field theory correspondence [112] lead to good agreements of m E− /m M+ for 1.5 < T /T pc < 3, as shown in Fig. 19 [28]. Further study is needed to clarify the meanings and implications of these results.

Summary
The WHOT-QCD Collaboration is pushing forward a series of projects to clarify the phase structure and thermodynamic properties of QCD matter at finite temperatures and densities, mainly using improved Wilson quarks. Wilson-type quarks do not have the theoretical lack of clarity of staggeredtype quarks, but require more computational resources. Thus, development and application of various computational techniques are mandatory. We have developed a T -integration method to calculate the equation of state in the fixed-scale approach, a Gaussian method based on the cumulant expansion to tame the sign problem, and a histogram method combined with the reweighting technique to explore the phase structure of QCD. Using them, we have made a series of studies in QCD at finite temperatures and densities with improved Wilson quarks. In particular, we have carried out the first study of finite-density QCD with two flavors of Wilson quarks and the first calculation of the equation of state with 2 + 1 flavors of Wilson quarks. We are currently extending the studies towards our final objective of 2 + 1 flavor QCD directly at the physical point.