The phase diagram of dense QCD
Abstract
Current status of theoretical researches on the QCD phase diagram at finite temperature and baryon chemical potential is reviewed with special emphasis on the origin of various phases and their symmetry breaking patterns. Topics include; quark deconfinement, chiral symmetry restoration, order of the phase transitions, QCD critical point(s), colour superconductivity, various inhomogeneous states and implications from QCDlike theories.
YITP1028, TKYNT1006
Contents
 1 Introduction
 2 QCD phase structure
 3 Order parameters for the QCD phase transition
 4 Chiral phase transition at finite temperature
 5 Chiral phase transition at finite baryon density
 6 Formation of the diquark condensate
 7 Inhomogeneous states
 8 Suggestions from QCDlike theories
 9 Summary and concluding remarks
1 Introduction
One of the most crucial properties in the nonAbelian gauge theory of quarks and gluons, Quantum Chromodynamics (QCD) [1], is the asymptotic freedom [2, 3]; the coupling constant runs towards a smaller value with increasing energy scale. It is hence a natural anticipation that QCD matter at high energy densities undergoes a phase transition from a state with confined hadrons into a new state of matter with onshell (real) quarks and gluons.
There are two important external parameters for QCD in equilibrium, the temperature and the baryon number density . (In the grand canonical ensemble, the quark chemical potential may be introduced as a conjugate variable to the quark number density .) Since the intrinsic scale of QCD is , it would be conceivable that the QCD phase transition should take place around or . Experimentally, the heavyion collisions (HIC) in laboratories provide us with a chance to create hot and/or denser QCD matter and elucidate its properties. In particular, the Relativistic HeavyIon Collider (RHIC) at Brookhaven National Laboratory (BNL) has conducted experiments to create hot QCD matter (a quarkgluon plasma or QGP in short) by the AuAu collisions with the highest collision energy [4]. The Large Hadron Collider (LHC) at CERN will continue experiments along the same line with higher energies [4, 5]. Exploration of a wider range of the QCD phase diagram with up to several times of the normal nuclear matter density may be carried out by lowenergy scan in HIC at RHIC as well as at the future facilities such as the Facility for Antiproton and Ion Research (FAIR) at GSI, the Nuclotronbased Ion Collider Facility (NICA) at JINR and the Japan Proton Accelerator Research Complex (JPARC) at JAERI.
In nature, the deep interior of compact stellar objects such as neutron stars would be the relevant place where dense QCD matter at low temperature is realized (see [6] for a review). In fact, continuous efforts have been and are being conducted in the observations of neutron stars to extract information of the equation of state of dense QCD. If the baryon density is asymptotically high, weak coupling QCD analyses indicate that the QCD ground state forms a condensation of quark Cooper pairs, namely the colour superconductivity (CSC). Since quarks have not only spin but also colour and flavour quantum numbers, the quark pairing pattern is much more intricate than the electron pairing in metallic superconductors.
In this review we will discuss selected topical developments in the QCD phase diagram with an emphasis on the phases at finite . For the topics not covered in the present article, readers may consult other reviews [7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19].
We organize this article as follows. In section 2, after a brief introduction of two key features of hot/dense QCD, i.e. the deconfinement and chiral restoration, we show a conjectured QCD phase diagram in – plane. In section 3 we introduce three major order parameters to characterize QCD matter at finite and , i.e. the Polyakov loop, the chiral condensate and the diquark condensate. In section 4 we summarize the current status on the chiral phase transition at finite with . In section 5 we go into more details about the phase transitions in the density region up to a moderate comparable to , which are classified according to different spontaneous chiral symmetry breaking (SB) patterns. We then address the situation with much greater than in section 6. Even though the theoretical understanding has not been fully settled down yet, section 7 is devoted to a review over conjectured inhomogeneous states of QCD matter as well as some alternative scenarios. In section 8 we look over suggestive results and implications from QCDlike theories. Section 9 is devoted to the summary and concluding remarks.
2 QCD phase structure
We will run through the QCD phase transitions and associated phase structure here before detailed discussions in subsequent sections.
2.1 Deconfinement and chiral restoration
A first prototype of the QCD phase diagram in  plane was conjectured in [20]. It was elucidated that one could give an interpretation of the Hagedorn’s limiting temperature in the Statistical Bootstrap Model (SBM) [21] as a critical temperature associated with a secondorder phase transition into a new state of matter. The weakly interacting quark matter at large due to the asymptotic freedom had been also recognized [22]. Historical summary of QCD phase diagram and its exploration in HIC experiments are given in the reviews [23, 24].

Deconfinement — In an early picture of hadron resonance gas at finite temperature [25], the density of (mostly mesonic) states as a function of the resonance mass is proportional to where is known from the Regge slope parameter. Such an exponentially growing behaviour of the density of states should be balanced by the Boltzmann factor in the partition function. When , the integration over becomes singular, so that plays a role of the limiting temperature (Hagedorn temperature) above which the hadronic description breaks down. This argument is applied to estimate the critical value of as well. The density of baryonic states, , is balanced by the Boltzmann factor , leading to the limiting temperature . We see equivalently that the critical at is given by .
If one bears in mind a simple bagmodel picture that hadrons are finitesize objects in which valence quarks are confined, it is conceivable to imagine that hadrons overlap with each other and start to percolate at the Hagedorn temperature [26, 27], which is an intuitive portrayal of quark deconfinement. Although the Hagedorn/percolation picture is useful for practical objectivization, it is necessary to develop a fieldtheoretical definition of the quarkdeconfinement in QCD. Global centre symmetry of pure gluonic sector of QCD gives such a definition as elucidated in section 3.1.

Chiral restoration — The QCD vacuum should be regarded as a medium with full of quantum fluctuations that are responsible for the generation of nonperturbative quark mass. In hot and dense energetic matter, quarks turn bare due to asymptotic freedom. Therefore, one may expect a phase transition from a state with heavy constituent quarks to another state with light current quarks. Such a transition is called chiral phase transition named after the underlying chiral symmetry of QCD. The QCD phase diagram at finite and was also conjectured from the point of view of chiral symmetry [28]. In this case, the order parameter is the chiral condensate which takes a value about in the vacuum and sets a natural scale for the critical temperature of chiral restoration. In the chiral perturbation theory (PT) the chiral condensate for two massless quark flavours at low temperature is known to behave as with the pion decay constant [29]. Although the validity of PT is limited to low temperature, this is a clear evidence of the melting of chiral condensate at finite temperature. At low baryon density, likewise, the chiral condensate decreases as [30, 31, 32] where is the  sigma term. (For higherorder corrections, see [33, 34].)
The chiral transition is a notion independent of the deconfinement transition. In section 3.2 we classify the chiral transition according to the SB pattern.
2.2 Conjectured QCD phase diagram
Figure 1 summarizes our stateoftheart understanding on the phase structure of QCD matter including conjectures which are not fully established. At present, relatively firm statements can be made only in limited cases – phase structure at finite with small baryon density () and that at asymptotically high density (). Below we will take a closer look at figure 1 from a smaller to larger value of in order.
Hadronquark phase transition at :
The QCD phase transition at finite temperature with zero chemical potential has been studied extensively in the numerical simulation on the lattice. Results depend on the number of colours and flavours as expected from the analysis of effective theories on the basis of the renormalization group together with the universality [35, 36]. A firstorder deconfinement transition for and has been established from the finite size scaling analysis on the lattice [37], and the critical temperature is found to be . For light flavours it is appropriate to address more on the chiral phase transition. Recent analyses on the basis of the staggered fermion and Wilson fermion indicate a crossover from the hadronic phase to the quarkgluon plasma for realistic , and quark masses [38, 39]. The pseudocritical temperature , which characterizes the crossover location, is likely to be within the range as summarized in section 4.2.
Even for the temperature above the system may be strongly correlated and show nonperturbative phenomena such as the existence of hadronic modes or preformed hadrons in the quarkgluon plasma at [28, 40] as well as at [41, 42, 43]. Similar phenomena can be seen in other strong coupling systems such as the hightemperature superconductivity and in the BEC regime of ultracold fermionic atoms [44].
QCD critical points:
In the density region beyond there is no reliable information from the firstprinciple lattice QCD calculation. Investigation using effective models is a pragmatic alternative then. Most of the chiral models suggest that there is a QCD critical point located at and the chiral transition becomes firstorder (crossover) for () for realistic , and quark masses [45, 46, 47, 48] (see the point E in figure 2). The criticality implies enhanced fluctuations, so that the search for the QCD critical point is of great experimental interest [49, 50].
There is also a possibility that the firstorder phase boundary ends at another critical point in the lower and higher region whose location we shall denote by as shown by the point F in figure 2. As discussed in section 6, the cold dense QCD matter with three degenerate flavours may have no clear border between superfluid nuclear matter and superconducting quark matter, which is called the quarkhadron continuity.
In reality, the fate of the above critical points (E and F) depend strongly on the relative magnitude of the strange quark mass and the typical values of and at the phase boundary.
Liquidgas phase transition of nuclear matter:
Since the nucleon mass is and the binding energy in isospinsymmetric nuclear matter is around , a nonvanishing baryon density of nuclear matter starts arising at at . At the threshold , the density varies from zero to the normal nuclear density . For the nuclear matter is fragmented into droplets with , so that is achieved on spatial average. This is a typical firstorder phase transition of the liquidgas type. The firstorder transition weakens as grows and eventually ends up with a secondorder critical point at as indicated by the point G in figure 2. Low energy HIC experiments indicate that and [51].
Quakyonic matter:
The Statistical Model is successful to reproduce the experimentally observable particle abundances at various and thus various . This model assumes a thermally equilibrated gas of noninteracting mesons, baryons and resonances for a given and [52, 53, 54]. Within the model description one can extract and from the HIC data by fitting the particle ratios. The accumulation of extracted points makes a curve on the – plane, which is called the chemical freezeout line.
The freezeout line is not necessarily associated with any of QCD phase boundaries. Nevertheless, there is an argument to claim that the sudden freezeout of chemical compositions should take place close to the phase transition [55]. Along the freezeout line the thermal degrees of freedom are dominated by mesons for . The higher becomes, the more baryons are excited. This indicates that there must be a transitional change at , where the importance of baryons in thermodynamics surpasses that of mesons. This happens around and according to the Statistical Model analysis.
It is interesting that such a phase structure is suggested from the large limit of QCD. When is large, quark loops are suppressed by as compared to gluon contributions [56, 57]. A finite baryon number density arises and the pressure grows of order once becomes greater than the lowest baryon mass . Such cold dense matter in the world is named quarkyonic matter [58]. Then, the phase diagram of large QCD consists of three regions separated by firstorder phase transitions, i.e. the confined, deconfined and quarkyonic phases. The meeting point of the three firstorder phase boundaries is the triple point whose remnant for finite is indicated by the point H in figure 2, as suggested in [59]. We will revisit the idea of quarkyonic matter in section 8.1.
Colour superconductivity:
If is asymptotically large, i.e. , the ground state of QCD matter can be analyzed in terms of the weakcoupling methods in QCD. Also we can count on the knowledge from condensed matter physics with quarks substituting for electrons. In this analogue between electrons in metal and quarks in quark matter, one may well anticipate that the ground state of QCD matter at low should form Cooper pairs leading to colour superconductivity (CSC) [13, 15, 18, 19, 60, 61]. Theoretical characterization will be further elucidated in section 3.3 and section 6.
There are many patterns of Cooper pairing and thus many different CSC
states. The search for the most stable CSC state still remains
unsettled except for or . It is in fact the
strange quark mass that makes the problem cumbersome. In the
intermediate density region particularly, the Fermi surface mismatch
of different quark flavours is given by
. When the gap energy is
comparable with , an inhomogeneous diquark condensation
may have a chance to develop energetically. Such an inhomogeneous CSC
state gives rise to a crystal structure with respect to ,
that is the crystalline CSC phase [62]. There are, at
the same time, various candidates over the crystalline CSC phase in
this intermediate density region and the true ground state has not
been fully revealed there (see section 7).
3 Order parameters for the QCD phase transition
QCD has (at least) three order parameters for quark deconfinement, chiral symmetry restoration and colour superconductivity as discussed below.
3.1 Polyakov loop and quark deconfinement
The Polyakov loop which characterizes the deconfinement transition in Euclidean spacetime is defined as [63, 64]
(1) 
which is an matrix in colour space. Here is the inverse temperature , and represents the path ordering. We will use to represent the traced Polyakov loop,
(2) 
Let us consider the centre of the colour gauge group : elements of the centre commute with all elements and can be written as with () and being an unit matrix. Under a nonperiodic gauge transformation of the following form; , the gauge fields receive a constant shift,
(3) 
so that the traced Polyakov loop transforms as . Because still keeps the periodicity in , such a nonperiodic gauge transformation still forms a symmetry of the gauge action. This is called centre symmetry [8, 35]. The quark action (with the quark field denoted by ) explicitly breaks centre symmetry because the transformed field does not respect the antiperiodic boundary condition any longer. Thus, centre symmetry is an exact symmetry only in the pure gluonic theory where dynamical quarks are absent or quark masses are infinitely heavy ().
The expectation value of the Polyakov loop and its correlation in the pure gluonic theory can be written as [65, 66, 67]
(4)  
(5) 
Here, the constant () independent of is an excess free energy for a static quark (antiquark) in a hot gluon medium. ^{1}^{1}1Strictly speaking, the expectation value in the pure gluonic theory should be defined by taking the limit after taking the thermodynamic limit , so that it takes a real value. Alternatively, one may use to define the free energy of a single quark. Also, is an excess free energy for an antiquark at and a quark at . ^{2}^{2}2Even when there are dynamical quarks, one may use the Polyakov loops, (4) and (5), to define the heavyquark free energies.
In the confining phase of the pure gluonic theory, the free energy of a single quark diverges () and the potential between a quark and an antiquark increases linearly at long distance ( with ), which leads to and . On the other hand, in the deconfined phase, the free energy of a single quark is finite (). Also the potential between a quark and an antiquark is of the Yukawa type at long distance with a magnetic screening mass [68, 69, 70],
(6) 
where is a dimensionless constant. Note that the glueball exchange with mass of which is smaller than the electric screening scale dominates over the longrange correlation at weak coupling. Therefore, and .
Confined (Disordered) Phase  Deconfined (Ordered) Phase  

Free Energy  
Polyakov Loop  
The qualitative behaviour of and is summarized in table 1. We see that can nicely characterize the state of matter in such a way that behaves like a magnetization in 3D classical spin systems [35, 8, 71]. The pure gluonic theory for and have been studied in lattice gauge simulations with the finitesize scaling analysis [72, 37]. It was shown that there is a secondorder phase transition for and a firstorder phase transition for : For the critical exponent is found to agree with the Ising model in accordance with the universality. For the GinzburgLandau free energy for the Polyakov loop has a cubic invariant, so that the firstorder transition can natually be induced [8, 37, 73]. Recent studies of the pure gluonic theory with , , , indicate that the transition is of first order for and becomes stronger as increases [74, 75]. This behavior is similar to that of the 3D state Potts model [76].
The idea to construct the deconfinement order parameter can be extended to a more general setup. Suppose that there is an arbitrary operator written in terms of gauge fields that is not invariant under centre transformation. (If we take as the quark propagator straight up along the imaginarytime direction, the following argument shall result in the standard definition of the Polyakov loop.) The expectation value of in the pure gluonic theory can be decomposed by the triality projection [77]; with
(7) 
where is defined in (3). Then, with can be an order parameter sensitive to the spontaneous breaking of centre symmetry. A particular choice [78, 79] (socalled the dual condensate) has some practical advantages [80].
Centre symmetry discussed so far is rigorously defined only in the pure gluonic system as already mentioned: In the presence of dynamical quarks, centre symmetry is broken explicitly, so that and always take finite values. This is analogous to the spin system under external magnetic field, in which the magnetization is always nonvanishing [81]. Nevertheless, may be estimated as with a hadronic mass scale () at low and thus the low phase approximately preserves centre symmetry.
3.2 Chiral condensate and dynamical breaking of chiral symmetry
In the QCD vacuum at , chiral symmetry is spontaneously broken, which is the source of hadron masses. It is a common wisdom that the chiral symmetry breaking is driven by the expectation value of an operator which transforms as under chiral symmetry. A simplest choice of the order parameter for the chiral symmetry breaking is a bilinear form called the chiral condensate,
(8) 
where colour and flavour indices of the quark fields are to be summed. If the above chiral condensate is nonvanishing even after taking the limit of zero quark masses, chiral symmetry is spontaneously broken according to the pattern with
(9)  
(10) 
which leads to massless NambuGoldstone bosons for . ^{3}^{3}3Strictly speaking, the quotient groups and need to be introduced to avoid double counting of the discrete centre group operation. A simplest example of this kind is [82]. The symmetry in the classical level of the QCD Lagrangian is broken down explicitly to in the quantum level. Then, the current is no longer conserved ( anomaly); ^{4}^{4}4The anomaly relation does not necessarily guarantee the absence of the massless NambuGoldstone boson. In the Schwinger model which is an exactly solvable QCD analogue, the problem is resolved by the KogutSusskind dipole ghosts.
(11) 
The righthand side of the above relation is nothing but the topological charge density. Thus, gauge configurations with nontrivial topology are microscopically responsible for the anomaly. In other words, the current could be approximately conserved if the gauge configurations are dominated by topologically trivial sectors. We will come back to this point later to address effective restoration of symmetry in the medium [83].
Note that (8) is not a unique choice for the order parameter [84, 85, 86, 87]. For example, one may consider the following fourquark condensate,
(12) 
If this is nonvanishing, the ground state breaks chiral symmetry to where corresponds to a discrete axial rotation. If the bilinear condensate (8) is nonzero, the fourquark condensate (12) takes a finite value in general. However, a nonzero value of (12) does not necessarily enforce a finite value of (8). In fact, if symmetry is left unbroken, (8) must vanish. As long as the Dirac determinant in QCD is positive definite, possibility of unbroken symmetry has been ruled out by the exact QCD inequality [85]. However, it is not necessary the case for finite with which the positivity of the Dirac determinant does not hold.
3.3 Diquark condensate and colour superconductivity
QCD at high baryon density shows a novel mechanism of spontaneous chiral symmetry breaking [88]. The fundamental degrees of freedom in the CSC phase with three colours and three flavours are the diquarks defined as
(13) 
where are flavour indices and are colour indices. Note that the charge conjugation matrix is necessary to make be Lorentz scalar. Then, is a triplet both in colour and flavour, so that it transforms in the same way as the quark field .
Under certain gauge fixing, one may consider the expectation values of these operators. They are called the diquark condensates as will be discussed further in section 6. Instead, we can also construct an analogue of (8) in terms of the diquarks,
(14) 
where colour and flavour indices are summed. This is a gaugeinvariant fourquark condensate which characterizes the spontaneous chiral symmetry breaking in CSC [13]. Unlike the bilinear operator in (8), the fourquark operator here keeps additional invariance corresponding to the reflections; and , which is broken down to due to the anomaly.
In the CSC phase, not only chiral symmetry but also symmetry associated with the baryon number conservation may be spontaneously broken. A colour singlet order parameter to detect such symmetry breaking can be
(15) 
The sixquark operator here breaks with its subgroup maintained. If this condensate is nonzero, there appears exactly massless NambuGoldstone boson because the baryon number symmetry is an exact symmetry in the QCD Lagrangian.
4 Chiral phase transition at finite temperature
The chiral phase transition at finite with has been and is being extensively studied by the renormalization group method near the critical point à la GinzburgLandauWilson and by the latticeQCD simulations. In this section we will briefly summarize the current status of these studies. (See [39] for further details.)
4.1 GinzburgLandauWilson analysis
If the phase transition is of second order or of weak first order, one may write down the freeenergy functional in terms of the order parameter field as a power series of . The large fluctuation of near the critical point is then taken into account by the renormalization group method. This is called the GinzburgLandauWilson approach. For chiral phase transition in QCD, the relevant order parameter field is a matrix in flavour space, . Under the flavour chiral rotation , transforms as . Then the GinzburgLandau free energy in three spatial dimensions () with full symmetry up to the quartic order in becomes [36, 47];
(16) 
Effects of temperature enter through the parameters , and . Note that is bounded from below as long as and are satisfied. The renormalization group analysis of (16) on the basis of the leadingorder expansion leads to a conclusion that there is no stable IR fixed point for [36]. This implies that the thermal phase transition described by (16) is of the fluctuationinduced first order for two or more flavours.
In QCD, however, there is anomaly and the correct chiral symmetry is for massless quarks. The lowest dimensional operator which breaks symmetry explicitly while keeping the rest of chiral symmetry is the Kobayashi–Maskawa–’t Hooft (KMT) term [89, 90, 91, 92];
(17) 
The coefficient , which is dependent in general, dictates the strength of anomaly. In the instanton picture [91, 92], is proportional to the instanton density , which is perturbatively evaluated as [7]
(18) 
with being a typical instanton size.
For the KMT term becomes a cubic invariant in the order parameter. Hence, leads to the chiral phase transition of first order. For , on the other hand, the KMT term becomes a quadratic invariant. Also the chiral symmetry in this case is . Such an effective theory with symmetry has a WilsonFisher type IR fixed point as long as the coefficient of the quartic term of is positive. Therefore, if the chiral phase transition of massless QCD is of second order, its critical exponents would be the same as those in the 3D effective theory according to the notion of universality. In Table 2 we summarize the GinzburgLandauWilson analysis from the chiral effective theory [36].
symmetric ()  Fluctuationinduced 1st order  1st order 

broken ()  2nd order [ universality]  1st order 
In the real world, none of quark is exactly massless: For example, –, – and at the renormalization scale of [93]. Therefore, it is useful to draw a phase diagram by treating quark masses as external parameters. This is called the Columbia plot [94] as shown in figure 3 where the isospin degeneracy is assumed (). The firstorder chiral transition and the firstorder deconfinement transition at finite are indicated by the leftbottom region and the righttop region, respectively. The chiral and deconfinement critical lines, which separate the firstorder and crossover regions, belong to a universality class of the 3D Ising model except for special points at or [95].
If the chiral transition is of second order for massless case, the chiral critical line meets the axis at (tricritical point) and changes its universality to for [96]. The tricritical point at is a Gaussian fixed point of the 3D model (that is, the critical dimension is not but at the tricritical point), so that the critical exponents take the classical (meanfield) values [97], which is confirmed in numerical studies of the chiral model [98].
4.2 Lattice QCD simulations
Although the critical properties expected from the GinzburgLandauWilson analysis discussed above are expected to be universal, the quantities such as the critical temperature and the equation of state depend on the details of microscopic dynamics. In QCD, only a reliable method known for microscopic calculation is the latticeQCD simulation in which the functional integration is carried out on the spacetime lattice with a lattice spacing and the lattice volume by the method of importance sampling. In latticeQCD simulations there are at least two extrapolations required to obtain physical results; the extrapolation to the continuum limit () and the extrapolation to the thermodynamic limit (). Therefore, lattice results receive not only statistical errors due to the importance sampling but systematic errors due to the extrapolations also.
For nearly massless fermions in QCD, there is an extra complication to reconcile chiral symmetry and lattice discretization; the Wilson fermion and the staggered fermion have been the standard ways to define light quarks on the lattice, while the domainwall fermion and the overlap fermion recently proposed have more solid theoretical ground although the simulation costs are higher. For various applications of latticeQCD simulations to the system at finite and , see a recent review [39].
Here we mention only two points relevant to the discussions below: (i) The thermal transition for physical quark masses is likely to be crossover as indicated by a starsymbol in figure 3. This is based on the finitesize scaling analysis using staggered fermion [38]. Confirmation of this result by other fermion formalisms is necessary, however. (ii) The (pseudo)critical temperature with different types of fermions and with different lattice spacings are summarized in figure 4. In view of these data with error bars, we adopt a conservative estimate at present; –. ^{5}^{5}5Possible uncertainties in stem from discretization errors, conversion from the lattice unit to the physical unit, and the prescription of defining . For crossover transition may depend on which susceptibility (either chiral susceptibility or Polyakov loop susceptibility ) are used, and also depend on dependent normalization of the susceptibilities. It has been clarified recently that improvement of the staggered action with less tastesymmetry breaking favours smaller value of [107, 108].
5 Chiral phase transition at finite baryon density
Let us now introduce the baryon chemical potential as an extra axis to the Columbia plot. In the socalled “standard scenario” the firstorder region in the lower left corner of the Columbia plot is elongated with increasing as written in the left panel of figure 5. If the physical point at is in the crossover region, there arises a critical chemical potential so that the system shows a firstorder transition for . That is, we find a QCD critical point at on the QCD phase diagram in the – plane. In the right of figure 5, the socalled “exotic scenario” is sketched, in which the size of the firstorder region shrinks as increases. In this case, if the physical point at is in the crossover region, it stays crossover for finite , so that there arises no critical point (at least for small ) in the QCD phase diagram in the – plane. In general the critical surface in figure 5 can have more complicated structure, which may allow for several QCD critical points in the – plane.
5.1 Lattice QCD at low baryon density
In the presence of , the QCD partition function on the lattice is written as
(19) 
where is the matrixvalued gauge field defined on the links. The YangMills action is denoted by with , while the quark contribution is denoted by the determinant of with being the Euclidean Dirac operator. The quark mass is real and positive. Let be an eigenfunction of as with an eigenvalue . If then is an antiHermitean operator and thus should be pure imaginary. Besides, is an independent eigenfunction with an eigenvalue , because of Hermiticity . Thus, we always have complexconjugate pairs and the determinant becomes real and positive;
(20) 
In this way, at , the integrand in (19) is positive definite and the importance sampling method works fine to evaluate the functional integral.
For the Hermiticity is replaced by
(21) 
Therefore is no longer real unless is pure imaginary. This requires us to deal with a severe cancellation of positive and negative numbers to evaluate the functional integral. This is the notorious sign problem whose difficulty grows exponentially as the lattice volume increases.
In the following we shall briefly look over some of approaches in the latticeQCD simulations at finite . For more details, see the reviews [109, 110, 111] and references therein.

Multiparameter reweighting method — An expectation value of an operator at finite can be formally rewritten in terms of an ensemble average at ;
(22) where is called the reweighting factor [112]. The gauge configurations generated at only occasionally sample the region where and are large (overlap problem). Also, they take complex values (sign problem). Therefore, such a simulation works only when and are small.
To have a better overlap, the reweighting may be generalized towards not only the direction but also the direction [113]:
(23) where the multiparameter reweighting factor is defined as with . Here should be chosen to maximize the overlap. The physical temperature and are implicitly related through the lattice spacing . Then, the reweighting along the phase boundary in the – plane (as sketched in the right panel of figure 6) would have a better chance to probe into larger .

Taylor expansion method — The full evaluation of the reweighting factor is not an easy task even on the computer. If is small enough, the righthand side of (22) can be expanded in terms of [114, 115, 116],
(24) The coefficients are written in terms of the quark propagator and can be simulated at . Because the expansion is based on (22), the original overlap problem translates into the convergence problem in higher order terms; the Taylor expansion makes sense for within a radius of convergence dictated by the singularity closest to the origin in the complex plane.

Imaginary chemical potential method — If is pure imaginary (which is here denoted by ), (21) reduces to the Hermiticity, so that there is no sign problem [117, 118]. Then, it is possible to perform lattice simulations to find physical observables by the analytic continuation back to the real chemical potential,
(25) The applicability of this method is bounded by singularities or periodicity. For imaginary chemical potential, plays a role of an angle variable which has naïve periodicity by . However, the partition function is a function of and a change in by can be absorbed by the centre transformation (3). Therefore, the actual period of is (RobergeWeiss (RW) periodicity) [119]. In the deconfinement phase at high , especially, there is a firstorder phase transition at a half of the RW point; for . Thus the method of imaginary chemical potential works up to at best, which is also confirmed in model studies [120].

Canonical ensemble method — In the thermodynamic limit the canonical ensemble with fixed particle number is equivalent with the grand canonical ensemble with fixed chemical potential (and the mean value of is specified). To convert the grand canonical to the canonical description, the Fourier transform in terms of the imaginary chemical potential is necessary [121, 122, 123, 117, 124, 125],
(26) In this case, the difficulty of the sign problem is transferred to the integration with respect to . This canonical approach works fine as long as the volume is not large, for which centre symmetry is forced to be restored for that is a multiple of [77]. It is a highly delicate procedure to take the correct thermodynamic limit in which not but should be kept fixed. Especially, the order of the phase transition is sensitive to how the thermodynamic limit is approached [126].

Density of states method — There are several variants of the density of states method depending on the choice of a variable to rewrite the partition function. Here we take an example of a phase of the Dirac determinant [127, 128]. (The plaquette or energy is useful as well [129, 130, 131].) The density of states in this case reads
(27) where the expectation value is taken by the phase quenched simulation in which is replaced by to avoid the sign problem. Using one can express the expectation value of an operator as
(28) with and
(29) Both and are quantities calculable in the lattice simulation. The difficulty of the sign problem is now translated into the precise determination of the density of states; (21) implies that the phase quenched simulation at with two degenerate flavours is identical to the simulation at finite isospin chemical potential ;
(30) Then, the phase quenched expectation value goes through a qualitative change as increases; an exotic phase with pion condensation appears for [132]. For small , both latticeQCD simulation [128] and an analysis in the PT [133, 134, 135] show that behaves as Gaussian, while becomes Lorenzian for in the PT. In the latter case, extremely precise determination of and is crucially important. So far, the density of states method is applied for relatively small values of combined with the Taylor expansion of the fermion determinant [128, 136].

Complex Langevin method — The sign problem arises from the importance sampling to deal with the multidimensional functional integral. Therefore, it may not appear in different quantization schemes other than the functional integral. The stochastic quantization [137, 138, 139, 140, 141] formulated in terms of the Langevin dynamics with a fictitious time is a promising candidate for such an alternative. If we consider a scalar field theory defined by an action , the Langevin equation reads
(31) where is a fictitious time and a Gaussian noise; and . The expectation value is taken over the distribution and the equilibrated value is obtained in the limit. At finite the action and noise are complex. The complex Langevin dynamics can evade the sign problem and correctly describe the system at finite density for some simple models, while there are also cases where it may fall into a wrong answer [142, 143]. Whether this method is applicable to dense QCD or not is still an open question.
It would be an important milestone if the latticeQCD simulation can show the existence/nonexistence of the QCD critical point in – plane. In [113, 144] the location of the LeeYang zero in the complex latticecoupling () plane has been investigated using the multiparameter reweighting for flavour staggered fermions. If the thermal transition is of first order, the LeeYang zero nearest to the real axis in a finite lattice box is expected to approach to the real axis as the lattice volume increases, while there is no such tendency for the crossover transition. It was then concluded that the QCD critical point is located at for physical quark masses [144]. It has been argued, however, in [145] that the sign problem at finite baryon density can fool the LeeYang zero near the real axis. It is also suggested to study the behaviour of a set of LeeYang zeros in the complex plane to make a firm conclusion.
If we assume that the convergence of the Taylor expansion in terms of is dictated by the singularity at the critical point, the radius of convergence is an indicator of the location of the critical point. In [116] it was reported that the coefficients up to sixth order yield and for flavour staggered fermions with . In the same way, with flavour staggered fermions, the radius of convergence is discussed in [146], though the results are not conclusive. The idea of the radius of convergence has been also tested in a chiral effective model coupled with the Polyakov loop, in which the exact location of the critical point and the higherorder coefficients are calculable in the meanfield approximation [147]. The result from the model test suggests that the relation between the convergence radius of the expansion and the location of the critical point is not clear.
In [148] the canonical partition function is extracted using the density of states method with an assumption of the Gaussian distribution of the phase of the Dirac determinant. In this case the standard Sshape curve in the – plane is a signature of the firstorder transition. The location of the critical point is estimated to be and for flavour staggered fermions with . In [125] the canonical ensemble method is used to study the Sshape and it was found and for flavour clover fermions with .
As shown in figure 5 one may consider the critical surface and its behaviour near to check whether the critical point is located in small region. In particular, on the symmetric line (), one can define the critical mass as a function of and make a Taylor expansion,
(32) 
The sign of and determines whether the critical surface behaves like the standard scenario or the exotic scenario in figure 5. Simulations with flavours of staggered fermions on a coarse lattice () with the imaginary chemical potential method show and [149], which is consistent with the exotic scenario in figure 5. This does not, however, exclude a possibility that the critical point exists for large value of . Confirmation of this result with different fermion formulations with smaller lattice spacing is necessary to draw a firm conclusion [150].
In short summary, in any method, the validity of QCD simulations at finite baryon density is still limited in the region at present because of the sign problem. Although some of the latticeQCD simulations suggest the existence of the QCD critical point in – plane, the results are to be taken with care if it is predicted at large .
5.2 Effective symmetry restoration
For the firstorder transition at finite is driven by the trilinear chiral condensates from the breaking KMTtype interaction with the strength as discussed in section 4.1. In the instanton calculation is proportional to the instanton density (18) which is screened by the medium at high and as
(33) 
in the oneloop calculation with the instanton size whose empirical value is [151]. It is thus expected that the breaking interaction is exponentially suppressed as and/or grow larger. If one performs the integration (which is not IR divergent thanks to the exponential suppression), the suppression factor would be or for by the dimensional reason. It is a nontrivial question whether the instantoninduced interaction drops by either exponential or power function [152]. The Columbia plot with an inmedium reduction of breaking should have a smaller region for the firstorder phase transition. This has been confirmed at [153] in the Polyakovloop extended chiral models such as the PNJL model [154, 155] and the PQM model [156]. Consequently the location of the QCD critical point at finite would be very sensitive to the strength [153, 157, 158]. The negative coefficients and in (32) might be attributed to the effective restoration at larger [159]; it has been found that the chiral model analysis with an exponential ansätz, , leads to a scenario similar to the right panel of figure 5.
5.3 QCD critical point search
In this subsection we summarize physical consequences of the critical point under the assumption that it exists in the QCD phase diagram.
Figure 7 is a schematic illustration of the shape of the effective potentials in the crossover and firstorder transitions. If the critical point is approached along the firstorder phase boundary, the associated potential shape looks like that of the secondorder phase transition around a nonvanishing order parameter. The critical properties at the critical point can be mapped to the 3D Ising model with symmetry with the (reduced) temperature and the external field . On the – plane in the QCD phase diagram the mapped direction is tangential to the firstorder phase boundary because symmetry at the critical point is locally preserved along this direction. The determination of the direction, on the other hand, requires microscopic calculations.
In the vicinity of the critical point along the direction the correlation length , the chiral condensate and the chiral susceptibility have the following scaling;
(34) 
where , and are known critical exponents in the 3D Ising model. Because chiral symmetry is explicitly broken by , the chiral condensate at the critical point is nonvanishing and takes a finite . In the same way, along the direction, they scale as