Phase diagram and critical properties within an effective model of QCD: the Nambu-Jona-Lasinio model coupled to the Polyakov loop

We investigate the phase diagram of the so-called Polyakov--Nambu--Jona-Lasinio model at finite temperature and non-zero chemical potential with three quark flavors. Chiral and deconfinement phase transitions are discussed and the relevant order-like parameters are analyzed. The results are compared with simple thermodynamic expectations and lattice data. We present the phase diagram in the ($T,\,\mu_B)$ plane, paying special attention to the critical end point: as the strength of the flavor-mixing interaction becomes weaker, the critical end point moves to low temperatures and can even disappear.


I. INTRODUCTION
Symmetries play a fundamental role in physics. In fact, the modern fundamental physics is dominated by considerations underlying symmetries which can be exact, approximate (explicitly broken) and (spontaneously) broken. A special role is played by gauge or local symmetries which lead to the description of the "real world" as the so-called local/gauge theories with spontaneously broken symmetries.
The concept of spontaneous symmetry breaking has been transferred from condensed matter physics to quantum field theory by Nambu [1]. It has been introduced in particle physics on the grounds of an analogy with the breaking of (electromagnetic) gauge symmetry in the theory of superconductivity by Bardeen, Cooper and Schrieffer (the so-called BCS theory). The application of spontaneous symmetry breaking to particle physics in the 1960s and successive years led to profound physical consequences and played a fundamental role in the edification of several models of elementary particles.
Spontaneous breaking of chiral symmetry, in particular, is known to govern the low-energy properties of hadrons [2][3][4]. Some QCD-like models have been proposed before the advent of QCD, and the phenomenon of spontaneous breaking of chiral symmetry and Nambu-Goldstone theorem was established more than 40 years ago [1].
The Nambu-Jona-Lasinio (NJL) model was proposed in 1961 to explain the origin of the nucleon mass with the help of spontaneous breaking of chiral symmetry [5,6]. At that time, the model was formulated in terms of nucleons, pions and scalar sigma mesons. The introduction of the quark degrees of freedom and the description of hadrons by Eguchi and Kikkawa [7,8] in the chiral limit, where the bare quark mass is m 0 = 0, and a more realistic version with m 0 = 0 by Volkov and Ebert [9][10][11], initiate a very intensive activity by several research groups [12].
There is strong evidence that quantum chromodynamics (QCD) is the fundamental theory of strong interactions. Its basic constituents are quarks and gluons that are confined in hadronic matter. It is believed that at high temperatures or densities the hadronic matter should undergo a phase transition into a new state of matter, the quark-gluon plasma (QGP).
A challenge of theoretical studies based on QCD is to predict the equation of state, the critical point and the nature of the phase transition.
As the evolution of QCD at finite density/temperature is very complicated, QCD-like models, as for instance NJL type models, have been developed providing guidance and information relevant to observable experimental signs of deconfinement and QGP features.
In fact, there has been great progress in the understanding of the properties of matter under extreme conditions of density and/or temperature, where the restoration of symmetries (e.g., the chiral symmetry) and the phenomenon of deconfinement should occur. These extreme conditions might be achieved in ultrarelativistic heavy-ion collisions or in the interior of neutron stars. In this context, increasing attention has been devoted to the study of the modification of particles propagating in a hot or dense medium [13,14]. The possible survival of bound states in the deconfined phase of QCD [15][16][17][18][19][20][21][22] has also opened interesting scenarios for the identification of the relevant degrees of freedom in the vicinity of the phase transition [23][24][25]. Besides lattice calculations [26][27][28][29], high temperature properties of QCD can be studied, starting from the QCD Lagrangian, within different theoretical schemes, like the dimensional reduction [30,31] or the hard thermal loop approximation [32][33][34]. Actually both the above approaches rely on a separation of momentum scales which, strictly speaking, holds only in the weak coupling regime g ≪ 1. Hence they cannot tell us anything about what happens in the vicinity of the phase transition. On the other hand, a system close to a phase transition is characterized by large correlation lengths (infinite in the case of a second order phase transition). Its behavior is mainly driven by the symmetries of the Lagrangian, rather than by the details of the microscopic interactions.
Confinement and chiral symmetry breaking are two of the most important features of QCD. As already referred, chiral models like the NJL model [5,6,35,36] have been successful in explaining the dynamics of spontaneous breaking of chiral symmetry and its restoration at high temperatures and densities/chemical potentials. Recently, this and other type of models, together with an intense experimental activity, are underway to construct the phase diagram of QCD.
The actual NJL type models describe interactions between constituent quarks, giving the correct chiral properties, and offers a simple and practical illustration of the basic mechanisms that drive the spontaneous breaking of chiral symmetry, a key feature of QCD in its low temperature and density phase. In order to take into account features of both chiral symmetry breaking and deconfinement, static degrees of freedom are introduced in this Lagrangian through an effective gluon potential in terms of the Polyakov loop [37][38][39][40][41][42][43][44]. The coupling of the quarks to the Polyakov loop leads to the reduction of the weight of the quark degrees of freedom at low temperature as a consequence of the restoration of the Z Nc symmetry associated with the color confinement.
In first approximation, the behavior of a system ruled by QCD is governed by the symmetry properties of the Lagrangian, namely the (approximate) global symmetry SU L (N f ) × SU R (N f ), which is spontaneously broken to SU V (N f ), and the (exact) SU c (N c ) local color symmetry. Indeed, in NJL type models the mass of a constituent quark is directly related to the chiral condensate, which is the order parameter of the chiral phase transition and, hence, is non-vanishing at zero temperature and density. Here the system lives in the phase of spontaneously broken chiral symmetry: the strong interaction, by polarizing the vacuum and turning it into a condensate of quark-antiquark pairs, transforms an initially point-like quark with its small bare mass m 0 into a massive quasiparticle with a finite size. Despite their widespread use, NJL models suffer a major shortcoming: the reduction to global (rather than local) color symmetry prevents quark confinement.
On the other hand, in a non-abelian pure gauge theory, the Polyakov loop serves as an order parameter for the transition from the low temperature, Z Nc symmetric, confined phase (the active degrees of freedom being color-singlet states, the glueballs), to the high temperature, deconfined phase (the active degrees of freedom being colored gluons), characterized by the spontaneous breaking of the Z Nc (center of SU c (N c )) symmetry. With the introduction of dynamical quarks, this symmetry breaking pattern is no longer exact: nevertheless it is still possible to distinguish a hadronic (confined) phase from a QGP (deconfined) one.
In the PNJL model, quarks are coupled simultaneously to the chiral condensate and to the Polyakov loop: the model includes features of both chiral and Z Nc symmetry breaking. The model has proven to be successful in reproducing lattice data concerning QCD thermodynamics [43]. The coupling to the Polyakov loop, resulting in a suppression of the unwanted quark contributions to the thermodynamics below the critical temperature, plays a fundamental role for the analysis of the critical behavior.
One of the important features of the QCD phase diagram is the existence of a phase boundary in the (T, µ B ) that separates the chirally broken hadronic phase from the chirally symmetric QGP phase. Arguments based on effective model calculations suggest that the QCD phase diagram can exhibit a tricritical point (TCP)/critical end point (CEP) where the line of first order matches that of second order/analytical crossover [45][46][47][48][49][50]. The discussion about the existence and location of such critical points of QCD is a very important topic nowadays [51]. This paper is organized as follows. In Section II we analyze the main features and symmetries of QCD. In Section III we present the model and formalism starting with the deduction of the self consistent equations. We also obtain the equations of state and the response functions. The regularization procedure used in the model calculations is also included. Section IV is devoted to the study of the equation of state at finite temperature.
In Section V we study the phase diagram and the location of the critical end point. In Section VI we discuss the important role of the choice of the model parameters for the correct description of isentropic trajectories. In Section VII we analyze the effects of strangeness and anomaly strength on the location of the critical end point. In Section VIII we proceed to study the size of the critical region around the critical end point and its consequences for the susceptibilities. Finally, concluding remarks are presented in Section IX.

A. Quantum Chromodynamics
The Lagrangian of QCD is written as [52,53] where q is the quark field with six flavors (u, d, s, c, b, t), three colors (N c = 3) andm being the corresponding current quark mass matrix in flavor space (m = diag f (m u , m d , . . . )). The covariant derivative, incorporates the color gauge field A a µ (a = 1, 2, ..., 8) and is the gluon field strength tensor; t a is the Gell-Mann color matrix in SU (3) ([t a , t b ] = if abc t c , tr(t a t b ) = δ ab 2 ) and f abc are the corresponding antisymmetric structure constants. Finally, g is the QCD coupling constant.
The QCD Lagrangian is by construction symmetric under SU(3) gauge transformations in color space and because of the non-Abelian character of the gauge group QCD has some main features: it is a renormalizable quantum field theory [54] with a single coupling constant for both, the quark-gluon interactions and the gluonic self-couplings involving vertices with three and four gluons; it has confinement, i.e., objects carrying color like quarks and gluons do not exist as physical degrees of freedom in the vacuum. In addition, QCD is a theory that has asymptotic freedom [55,56], i.e., for large momenta, Q, or wavelengths of the order of 10 −1 fm (ultraviolet region) the couplings are weak and the quarks and gluons propagate almost freely. For low momenta, or wavelengths of about 1 fm (infrared region), it happens the opposite situation and the couplings are quite strong. The system is now highly non perturbative. According to this property, the attraction between two quarks grows indefinitely as they move away from each other. This implies that the interaction between quarks and gluons can not be treated perturbatively, making the perturbative treatment of QCD not applicable to describe hadrons with masses below ∼2 GeV.
The low energy regime is specially interesting once it is relevant to the study of hadronic properties as for example in low energy QCD and nuclear physics.
The non-perturbative structure of the vacuum is characterized by the existence of quark condensates, i.e., it is expected non-zero values for the scalar density qq , by the appearance of light pseudoscalar particles, who are identified with (quasi) Goldstone bosons [57,58], and also by the existence of gluon pairs [59].
On the other hand, if QCD has a mechanism that includes the confinement of quarks the mass parameters m i are not observable quantities. However, they can be estimated in terms of the masses of some hadronic observables through current algebra's methods. These masses are denominated by current quark masses to distinguish them from the constituent quark masses, which are effective masses generated by the spontaneous breaking of chiral symmetry in phenomenological models of quarks.

B. Chiral Symmetry Breaking
Due to the relevant role that spontaneous breaking of chiral symmetry plays in hadronic physics at low energies, this symmetry is one of most important symmetries of QCD. Here we will concentrate in the N f = 3 case. In the chiral limit, i.e. m u = m d = m s = 0, QCD is chiral invariant which means that the QCD Lagrangian (1)   of symmetries: These symmetries are presented in Table I where it is also possible to see the transformations under which the Lagrangian is invariant, the currents that are conserved according to Noether's theorem and the respective manifestations of the symmetries in nature.
The SU V (3) and U V (1) symmetries ensure the conservation of isospin and baryon number, respectively, while the SU A (3) and U A (1) symmetries are transformations that involve the γ 5 matrix and therefore alter the parity of the state in which they operate. For the sake of uniformity throughout the text we will designate by chiral symmetry the SU A (3) symmetry and by axial symmetry the U A (1) symmetry.
From the experimental point of view the manifestation of chiral symmetry would be the existence of parity doublets, i.e., a multiplet of particles with the same mass and opposite parity for each multiplet of isospin (the chiral partners), in the hadronic spectrum; that situation is not verified. Similarly, if the symmetry U A (1) would manifest itself, the existence of a partner with opposite parity to each hadron should be observed experimentally. As neither of these situations is observed in the hadron spectrum, these symmetries must be somehow broken.
Concerning the SU A (3) symmetry, the theory must contain a mechanism for the spontaneous breaking of chiral symmetry, which represents a transition to an asymmetric phase. This is closely related to the existence of non-zero quark condensates, qq , which are not invariant under SU A (3) transformations and therefore act as order parameters for the spontaneously broken chiral symmetry.
According to the Goldstone theorem, the spontaneous breaking of a continuous global symmetry implies the existence of a particle with zero mass, the Goldstone boson. In the case we are considering, the symmetry breaking is closely related to the appearance of eight degenerate Goldstone bosons with zero mass. As matter of fact the pions where the first mesons associated to Goldstone bosons due to their small mass. Indeed, if compared to the mass of the nucleon one has M π /M N = 0.15. To reproduce the meson spectrum it is also necessary that the theory incorporates a mechanism which explicitly breaks the chiral symmetry: the Lagrangian must include perturbative terms that break ab initio the symmetry allowing the lifting of the degeneracy in the pseudoscalar mesons spectrum. These  Now we will analyze the U A (1) symmetry. It is known that in the chiral limit, and at the classical level, L QCD is also invariant under the axial transformation (see Table I).
As already mentioned, the absence in nature of chiral partners with the same mass, related with the U A (1) symmetry, opens the possibility that U A (1) symmetry is also spontaneously broken, similarly to what happens with the SU A (3) symmetry. Consequently, there should exist another pseudoscalar Goldstone boson. S. Weinberg estimated the mass of this particle, outside the chiral limit, at about √ 3M π [60].
This interaction "absorbs" N f left helicity fermions and converts them to right handed ones (and conversely). In the following, we will take θ inst = 0.
In this context, as suggested by 't Hooft, the instantons can play a crucial role in breaking explicitly the U A (1) symmetry giving to the η ′ a mass of about 1 GeV outside the chiral limit. This implies that the mass of η ′ has a different origin than the other masses of the pseudoscalar mesons and can not be seen as the missing Goldstone boson due to spontaneous breaking of U A (1) symmetry. Consequently, the U A (1) anomaly is very important once it is responsible for the flavor mixing effect that removes the degeneracy among several mesons. In this Section, following the arguments given in [38,39], we discuss how the deconfinement phase transition in a pure SU(N c ) gauge theory can be conveniently described through the introduction of an effective potential for the complex Polyakov loop field, which we define in the following.
Since we want to study the SU(N c ) phase structure, first of all an appropriate order parameter has to be defined. For this purpose the Polyakov line is introduced. In the above, A 4 = iA 0 is the temporal component of the Euclidean gauge field ( A, A 4 ), in which the strong coupling constant g has been absorbed, P denotes path ordering and the usual notation β = 1/T has been introduced with the Boltzmann constant set to one (k B = 1).
The Polyakov line L ( x) can be described as an operator of parallel transport of the gauge field A 4 ( x, τ ) into the direction τ . One way to understand why this quantity is indeed a parameter that can distinguish between confined or deconfined phase is to consider the two extremal behavior of L: a color field A 4 at the point ( x, τ ) will be transformed into the color field L × A 4 after being transported into the direction τ . If L → 1 it means that nothing affects the propagation of the field: the medium is in its deconfined phase. At the contrary L → 0 indicates that the color field cannot propagate in the medium: it is confined.
Another way to see this point is to consider the variation of free energy when an infinitely massive (hence static) quark (that acts as a test color charge) is added to the system. To this purpose let us introduce the Polyakov loop. When the theory is regularized on the lattice, it reads: and it is a color singlet under SU(N c ), but transforms non-trivially, like a field of charge one, under Z Nc . Its thermal expectation value can then be chosen as an order parameter for the deconfinement phase transition [63][64][65]: in the usual physical interpretation [66,67], l( x) is related to the change of free energy occurring when a heavy color source in the fundamental representation is added to the system. One has: In the Z Nc symmetric phase, l( x) = 0, implying that an infinite amount of free energy is required to add an isolated heavy quark to the system: in this phase color is confined.
In the case of the SU(3) gauge theory, the Polyakov line L( x) gets replaced by its gauge covariant average over a finite region of space, denoted as L( x) [38,39]. Note that L( x) in general is not a SU(N c ) matrix. The Polyakov loop field, is then introduced.
Following the Landau-Ginzburg approach, a Z 3 symmetric effective potential is defined for the (complex) Φ field, which is conveniently chosen to reproduce, at the mean field level, results obtained in lattice calculations [38,39,43]. In this approximation one simply sets the Polyakov loop field Φ( x) equal to its expectation value Φ = const., which minimizes the potential.
Concerning the effective potential for the (complex) field Φ, different choices are available in the literature [68][69][70]; the one proposed in [69] (see Equation (10)) is known to give sensible results [69,71,72] and will be adopted in our parametrization of the PNJL model that will be presented in Section III. In particular, this potential reproduces, at the mean field level, results obtained in lattice calculations as it will be shown. The potential reads: where The effective potential exhibits the feature of a phase transition from color confinement (T < T 0 , the minimum of the effective potential being at Φ = 0) to color deconfinement (T > T 0 , the minima of the effective potential occurring at Φ = 0).
The parameters of the effective potential U are given in Table II. These parameters have been fixed in order to reproduce the lattice data for the expectation value of the Polyakov loop and QCD thermodynamics in the pure gauge sector [73,74].
The parameter T 0 is the critical temperature for the deconfinement phase transition within a pure gauge approach: it was fixed to 270 MeV, according to lattice findings. Different criteria for fixing T 0 may be found in the literature, like in [75], where an explicit N f dependence of T 0 is presented by using renormalization group arguments. However, it is It should be noticed that the NJL parameters and the Polyakov potential ones are not on the same footing. Whereas the NJL parameters are directly related in a one to one correspondence with a physical quantity, the Polyakov loop potential is there to insure that the pure gauge lattice expectations are recovered. Hence the potential for the loop can be viewed as a unique but functional parameter. The details of this function are not very important to study the thermodynamic as soon as the potential reproduces the lattice results. In order to clarify this point, we remember that in a previous calculations in the N f = 2 case [76] we used two kind of potentials and obtained very small differences in what concerns thermodynamics; however, if one calculates susceptibilities with respect to Φ the log potential used here has to be preferred in order not to have unphysical results. The only true parameter is the pure gauge critical temperature T 0 that fixes the temperature scale of the system. However we would like to point out that it is expected in the Landau-Ginzburg framework that a characteristic temperature for a phase transition will not be a prediction: one needs to fix the correct energy scale somehow and it is the role of this parameter. Hence we will allow ourselves to change this parameter in the calculation of some observables in order to compare our results with lattice QCD expectations.
Finally we want to stress that, contrarily to the full Landau-Ginzburg effective field approach, the Polyakov loop effective field is not a dynamical degree of freedom due to the the lack of dynamical term in the Polyakov loop potential. Hence it is a background gauge field in which quarks will propagate. Anyway the potential mimics a pressure term that has the correct magnitude and temperature behavior in order to get the Stefan-Boltzman limit for the gluonic degrees of freedom and that explains the success of the model. Near the critical temperature, T c , the energy density (and other thermodynamic quantities) show a strong growth, signaling the transition from a hadronic resonance gas to a matter of deconfined quarks and gluons. As matter of fact, the rapid rise of the energy density is usually interpreted to be due to deconfinement i.e., liberation of many new degrees of freedom. In the limit where quark masses are infinitely heavy the order parameter for the deconfinement phase transition is the Polyakov loop. A rapid change in this quantity is also an indication for deconfinement even in the presence of light quarks [77].
Concerning chiral symmetry, and considering the chiral limit, it is expected a chiral transition with the corresponding order parameter being the quark condensate: the quark condensate vanishes at the critical temperature T c and a genuine phase transition takes place. Even away from the chiral limit, when the quark masses are finite, it is expected, in the transition region, a "crossover" where the quark condensates rapidly drop indicating a partial restoration of the chiral symmetry [77].
This rises the interesting question whether the restoration of chiral symmetry and the transition to the QGP occurs at same time. In [78] it was proposed that the evidence of restoration of chiral symmetry is a sufficient condition to demonstrate the existence of a new state of matter, but it is not a necessary condition for the discovery of the QGP. In fact, as already pointed out, most of lattice results shows a tendency for the restoration of chiral symmetry to happen simultaneously with deconfinement, but this issue is not definitively determined from a theoretical point of view.
At finite temperature and chemical potential the most common three-flavor phase diagram shows a first order boundary of the chiral phase transition separating the hadronic and quark phases. This first-order line starts at non-zero chemical potential and zero temperature and will finish in a point, the critical end point (T CEP , µ CEP B ), where the phase transition is of second order. As the temperature increases and the chemical potential decreases, the QCD phase transition becomes a crossover.
Although recent lattice QCD results by de Forcrand and Philipsen question the existence of the CEP [79][80][81], this critical point of QCD, proposed at the end of the eighties [45][46][47][48], is still a very important subject of discussion nowadays [51]. The search of the QCD CEP is one of the main goals in "Super Proton Synchrotron" (SPS) at CERN [82] and in the next phase of the "Relativistic Heavy Ion Collider" (RHIC) running at BNL.
When matter at high density and low temperature is considered, it is expected that this type of matter is a color superconductor where pairs of quarks condense ("diquark conden-  [36,83].
Recently it was argued that some features of hadron production in relativistic nuclear collisions, mainly at SPS-CERN energies, can be explained by the existence of a new form of matter beyond the hadronic matter and the QGP: the so-called quarkyonic matter [70,[84][85][86]. It was also suggested that these different types of matter meet at a triple point in the QCD phase diagram: both the hadronic matter, the QGP, and the quarkyonic matter all coexist [87]. It should also be noted that atomic nuclei by themselves represent a system at finite density and zero temperature. At normal nuclear density it is estimated that the quark condensate undergoes a reduction of 30% [93] so that the resulting effects in the chiral order parameter can be measured in beams of hadrons, electrons and photons in nuclear targets.
If the evidence for deconfinement can be found experimentally, then the search for manifestations of chiral symmetry restoration will be, in the near future, one of the main goals once the investigation of the properties of matter can provide a clear evidence for changes in the fundamental vacuum of QCD, with far-reaching consequences [78].

III. THE PNJL MODEL WITH THREE FLAVORS
The Polyakov-Nambu-Jona-Lasinio (PNJL model) that we want to build is an effective model for QCD, written in term of quark degrees of freedom. In the next paragraphs we will present the ingredients we need to build our model namely: a massive Dirac Lagrangian together with a four quark chirally invariant interaction (the original NJL Lagrangian with a small mass term that breaks explicitly the chiral symmetry, has observed in the mass spectrum); the so-called 't Hooft interaction that reproduces the interaction of quarks with instantons and finally a Polyakov loop potential that mimic the effect of the pure gauge (Yang-Mills) sector of QCD on the quarks.

A. Nambu-Jona-Lasinio Model with Anomaly and Explicit Symmetry Breaking
Phase transitions are usually characterized by large correlation lengths, i.e., much larger than the average distance between the elementary degrees of freedom of the system. Effective field theories then turn out to be a useful tool to describe a system near a phase transition.
In particular, in the usual Landau-Ginzburg approach, the order parameter is viewed as a field variable and for the latter an effective potential is built, respecting the symmetries of the original Lagrangian. The existence of a phase transition between two sectors where the chiral symmetry is spontaneously broken or restored (transition associated to the quark condensate that acts as an order parameter) and the Ginzburg-Landau theory suggests to use the symmetry motivated NJL Lagrangian [35,36,59,94] for the description of the coupling between quarks and the chiral condensate in the scalar-pseudoscalar sector. The associated Lagrangian which complies with the underlying symmetries of QCD described in the previous section reads:  The current quark mass matrixm is in general non degenerate and explicitly breaks chiral symmetry SU L (3) × SU R (3) to SU f (3) or its subgroup. In the following we will take m u = m d hence keeping exact the isospin symmetry.
Let us notice that NJL has no build-in confinement and is non-renormalizable thus it requires the introduction of a cutoff parameter Λ. temporal background gauge field Φ [70,96,97]. The Lagrangian reads: The covariant derivative is defined as D µ = ∂ µ − iA µ , with A µ = δ µ 0 A 0 (Polyakov gauge); in Euclidean notation A 0 = −iA 4 . The strong coupling constant g is absorbed in the definition of A µ (x) = gA µ a (x) λa 2 , where A µ a is the (SU c (3)) gauge field and λ a are the (color) Gell-Mann matrices.
At T = 0, it can be shown that the minimization of the grand potential leads to Φ =Φ = 0. So, the quark sector decouples from the gauge one, and the model is fixed as referred in the previous subsection.
Some remarks are in order concerning the applicability of the PNJL model. It should be noticed that in this model, beyond the chiral point-like coupling between quarks, the gluon dynamics is reduced to a simple static background field representing the Polyakov loop (see details in [44,68]). This scenario is expected to work only within a limited range of temperatures, since at large temperatures transverse gluons are expected to start to be thermodynamically active degrees of freedom and they are not taken into account in the PNJL model. We can assume that the range of applicability of the model is roughly limited to T ≤ (2−3)T c , since, as concluded in [98], transverse gluons start to contribute significantly for T > 2.5 T c , where T c is the deconfinement temperature. We will work in the (Hartree) mean field approximation. In this context, the quarks can be seen as free particles whose bare current masses m i are replaced by the constituent (or dressed) masses M i .
The quark propagator in the constant background field A 4 is then: In the above, p 0 = iω n and ω n = (2n + 1)πT is the Matsubara frequency for a fermion.
Within the mean field approximation, it is straightforward (see Ref. [35]) to obtain effective quark masses from the Lagrangian (16); these masses are given by the so-called gap equations: where the quark condensates q i q i , with i, j, k = u, d, s (to be fixed in cyclic order), have to be determined in a self-consistent way as explained in the Appendix.
The PNJL grand canonical potential density in the SU f (3) sector can be written as where E i is the quasi-particle energy for the quark i: E i = p 2 + M 2 i , and z + Φ and z − Φ are the partition function densities.
The explicit expression of z + Φ and z − Φ are given by: where E of the loop to the 1-and 2-particles Boltzmann factor will be discussed in Section VIII.

D. Equations of State and Response Functions
The equations of state can be derived from the thermodynamical potential Ω(T, µ). This allows for the comparison of some of the results with observables that have become accessible in lattice QCD at non-zero chemical potential.
As usual, the pressure p is defined such as its value is zero in the vacuum state [36] and, since the system is uniform, we have where V is the volume of the system.
The relevant observables are the baryonic density and the (scaled) "pressure difference" given by Due to the relevance for the study of the thermodynamics of matter created in relativistic heavy-ion collisions, it is interesting to perform an analysis of the isentropic trajectories. The equation of state for the entropy density, s, is given by and the energy density, ǫ, comes from the following fundamental relation of thermodynamics The energy density and the pressure are defined such that their values are zero in the vacuum state [36].
The baryon number susceptibility, χ B , and the specific heat, C, are the response of the baryon number density, ρ B (T, µ), and the entropy density, s(T, µ), to an infinitesimal variation of the quark chemical potential µ and temperature, given respectively by: These second order derivatives of the pressure are relevant quantities to discuss phase transitions, mainly the second order ones.

E. Model Parameters and Regularization Procedure
As already seen, the pure NJL sector involves five parameters: the coupling constants g S and g D , the current quark masses m u = m d , m s and the cutoff Λ defined in Section III.
These parameters are determined in the vacuum by fitting the experimental values of several physical quantities. We notice that the coupling constant g S and g D and the parameter Λ are correlated. For instance, if we increase g S in order to provide a more significant attraction between quarks, we must also increase the cutoff Λ in order to insure a good agreement with experimental results. In addition, the value of the cutoff itself does have some impact as far as the medium effects in the limit T = 0 are concerned. Also the strength of the coupling g D has a relevant role in the location of the CEP, as it will be discussed.
In fact, the choice of the parametrization may give rise to different physical scenarios at T = 0 and µ B = µ = 0 [36], even if they give reasonable fits to hadronic vacuum observables and predict a first order phase transition. The set of parameters we used insures the stability conditions and, consequently, the compatibility with thermodynamic expectations.
On the other hand, the regularization procedure, as soon as the temperature effects are considered, has relevant consequences on the behavior of physical observables, namely on the chiral condensates and the meson masses [99]. Advantages and drawbacks of these regularization procedures have been discussed within the NJL [99] and PNJL [76,100] models. We remind that one of the drawbacks of the regularization that consists in putting a cutoff only on the divergent integrals is that, at high temperature, there is a too fast decrease of the quark masses that become lower than their current values. This leads to a non physical behavior of the quark condensates that, after vanishing at the temperature (T ef f ) where constituent and current quark masses coincide, change sign and acquire nonzero values again. To avoid this unphysical effects we use the approximation of imposing by hand the condition that, above T ef f , M i = m i and q i q i = 0.

A. Characteristic Temperatures
At zero temperature and chemical potential, the chiral symmetry of QCD is explicitly broken. It is expected that chiral symmetry will be restored at high temperature; hence a phase transition occurs separating the low and the high temperature regions. This phenomena may be realized in high energy heavy ion collision experiments.
In this case (µ B = 0 and T = 0), where a crossover transition occurs in the PNJL model, a unique critical temperature cannot be defined, but one of many characteristic temperatures for the transition may be used. Of course, these temperatures should coincide in the limit where the transition becomes of second order (in the chiral limit, for example, when one is concerned with the restoration of chiral symmetry).
We start our analysis by identifying the characteristic temperatures which separate the different thermodynamic phases in the PNJL model [44], using the regularization that allows high momentum quark states at finite temperature (the cutoff is used only in the integrals that are divergent and Λ → ∞ in the ones that are already convergent because of the thermal distributions). Let us analyze the behavior of the quark masses and the field Φ. The characteristic temperature related to the deconfinement phase transition is T Φ c and the chiral phase transition characteristic temperature, T χ c , signals partial restoration of chiral symmetry. These temperatures are chosen to be, respectively, the inflexion points of the "quasi" order parameter Φ and of the chiral condensate q u q u ; as in [68], we define T c as the average of the two transition temperatures, T Φ c and T χ c . As shown in [100], the present regularization lowers the characteristic temperatures and decreases the gap T Φ c − T χ c , leading therefore to better agreement with lattice results.
Let us remark that this section is also devoted to the study of thermodynamic quantities (pressure, energy per particle and entropy) for which a rescaling of the temperature T 0 is needed, in order to get a better agreement with lattice results for these quantities. Therefore, following the argumentation presented in [68], here we will use the reduced temperature T c by rescaling the parameter T 0 from 270 to 210 MeV (let us stress that this rescaling is only done for the remainder of this section). Results for the characteristic temperatures with T 0 = 210 MeV and T 0 = 270 MeV are shown in Table III.
It should be noticed that, when T 0 = 210 MeV is used, we loose the almost perfect coincidence of the chiral and deconfinement transitions (they are shifted relative to each other by about 32 MeV) and we have T c = 187 MeV within the range expected from lattice calculations [77]. However, the behavior of the relevant physical quantities is qualitatively the same whether T 0 = 270 MeV or T 0 = 210 MeV.
An interesting point to be noticed in Figure 2 to a faster decrease of the quark masses around T c and the present regularization enhances this effect, even at temperatures higher than T c . At T ef f = 345 MeV the regularization is responsible for the full restoration of the chiral symmetry that was dynamically broken: the quark masses go to their current values and the quark condensates vanish (see Figure   2, left and right panels). As already shown in the framework of the pure NJL model [99] and in the PNJL model [100], the effect of allowing high momentum quark states is stronger for the strange quark mass. Indeed, with the conventional regularization, the non-strange constituent quark mass at high temperature is already very close to its current mass, the new regularization only enhancing this behavior. Differently, the strange quark mass is always far from its current value, unless we allow high momentum quarks to be present and, in that case, its constituent mass decreases very substantially and comes to its current value.
As shown in [100], at T ef f the behavior of some given observables signals the effective restoration of chiral and axial symmetry: the masses of the meson partners of both chiral and axial symmetry are degenerated and the topological susceptibility vanishes as we can see from Figure 2 (lower right panel).

B. Thermodynamic Quantities
In the limit of vanishing quark chemical potential, significant information on the phase structure of QCD at high temperature is obtained from lattice calculations. The transition to the phase characteristic of this regime is related with chiral and deconfinement transitions which are the main features of our model calculation.
In Figure 3, we plot the scaled pressure, the energy and the entropy as functions of the temperature compared with recent lattice results (see Reference [77]). Since the transition to the high temperature phase is a rapid crossover rather than a phase transition, the pressure, the entropy and the energy densities are continuous functions of the temperature.
We observe a similar behavior in the three curves: a sharp increase in the vicinity of the transition temperature and then a tendency to saturate at the corresponding ideal gas limit.
Asymptotically, the QCD pressure for N f massless quarks and (N 2 c − 1) massless gluons is given (µ B = 0) by: where the first term denotes the gluonic contribution and the second term the fermionic one.
The results follow the expected tendency and go to the free gas values (or Stefan-Boltzmann limit), a feature that was also found with this type of regularization in the context of the SU(2) PNJL model [68,103,104]. The inclusion of the Polyakov loop effective potential U(Φ,Φ; T ) (it can be seen as an effective pressure term mimicking the gluonic degrees of freedom of QCD) is required to get the correct limit (indeed in the NJL model the ideal gas limit is far to be reached due to the lack of gluonic degrees of freedom).
The inclusion of the Polyakov loop and the regularization procedure are essential to obtain the required increase of extensive thermodynamic quantities, insuring the convergence to the Stefan-Boltzmann (SB) limit of QCD [105]. Some comments are in order concerning the role of the regularization procedure for T > T c . In this temperature range, due to the presence of high momentum quark states, the physical situation is dominated by the significant decrease of the constituent quark masses by the qq interactions. This allows for an ideal gas behavior of almost massless quarks with the correct number of degrees of freedom.
Let us notice that, just below T c , the pressure and the energy fail to reproduce the lattice points: for example there is a small underestimation of the pressure and energy in the model calculations. It is known that the lack of mesonic correlations in the PNJL model is responsible for, at least, a fraction of this discrepancy. As matter of fact, in [106]  is still a very important subject of discussion nowadays [51].
We remember that the TCP separates the second order transition at high temperature and low chemical potential, from the first order transition at high chemical potential and low temperature. If the second order transition is replaced by a smooth crossover, a CEP which separates the two lines is found. In order to determine and elucidate the nature of the phase transition the relevant thermodynamic quantities are studied, starting with zero temperature and finite chemical potential.
A. Phase Transition at Zero Temperature by the minimum of the thermodynamic potential. When stable and metastable solutions give the same value for the thermodynamic potential, the phase transition occurs as illustrated in Figure 4 (left panel). The phase of broken symmetry is realized for µ B < µ cr B and the "symmetric" phase is realized for µ B > µ cr B . At this crossing point of the curve, the two phases are in thermal and chemical equilibrium (Gibbs criteria). The baryon density, represented in Figure 4 (right panel) as function of µ B , is given by the slope of the curve −Ω as indicated by Equation (23). This quantity is plotted in Figure 4 (right panel) which shows that the condition of thermodynamic stability (∂(−Ω)/∂µ B > 0) is violated by the portion of the curve with negative curvature (unstable phase).
This information can be complemented by the behavior of the pressure/energy per particle as a function of the baryonic density. To this purpose let us analyze the curve at T = 0 in Figure 5. The pressure has three zeros that correspond to the extrema of the energy per particle. The third zero of the pressure, at ρ B = 2.36ρ 0 , corresponds to an absolute minimum of the energy (see Figure 5 right panel). This is an important point of the model calculation, and the set of parameters is chosen in order to insure such a condition. In fact, the occurrence of an absolute minimum of the energy allows for the existence of finite droplets in mechanical equilibrium with the vacuum at zero-pressure (P = 0). For densities above a critical value, ρ cr B = 2.36ρ 0 , the system returns to a uniform gas phase. The equilibrium configuration for densities 0 < ρ B < ρ cr B is, therefore, a mixed phase. In view of the behavior above described, we can conclude that, for T = 0, the uniform non-zero density phase will break up into stable droplets, with zero pressure and density ρ cr B = 2.36ρ 0 , in which chiral symmetry is partially restored, surrounded by a nontrivial vacuum with ρ B = P = 0 (see also [36,[107][108][109][110][111]  To determine the range of the crossover region we choose to define it as the interval between the two peaks around the zero of ∂ 2 q u q u /∂T 2 . This area is presented in gray in Figure 6 (left panel). We can see that as µ B increases the area where the crossover takes place is narrowed until we reaches the CEP.
Finally, we will focus again on the energy per baryon. In Figure 5 (right panel), we plot the density dependence of the energy per baryon at different temperatures. We observe that the two points (zero of the pressure and minimum of the energy density) are not the same at finite temperature. In fact, as can be seen from Figure 5  coincide with the minimum of the energy per particle (see Figure 5, right panel).
The arguments just presented allow to clarify the difference between confined quark matter (in hadrons) and bounded quark matter (droplets of quarks).
This pattern of phase transition is similar to the liquid-gas transition in nuclear matter, a consequence of the fact that nuclear matter assumes its ground state at a non-vanishing baryon density ρ B ∼ 0.17 fm −3 when T = 0.

VI. NERNST PRINCIPLE AND ISENTROPIC TRAJECTORIES
The isentropic lines contain important information about the conditions that are supposed to be realized in heavy ion collisions. Most of the studies on this topic have been done with lattice calculations for two flavor QCD at finite µ [113] but there are also studies using different type of models [112,114,115]. Some model calculations predict that in a region around the CEP the properties of matter are only slowly modified as the collision energy is changed, as a consequence of the attractor character of the CEP [116].
Our numerical results for the isentropic lines in the (T, µ B ) plane are shown in Figure 7.
We start the discussion by analyzing the behavior of the isentropic lines in the limit T → 0.
As already analyzed in Section V, our convenient choice of the model parameters allows a better description of the first order transition than other treatments of the NJL (PNJL) model. This choice is crucial to obtain important results: the criterion of stability of the quark droplets [36,109] is fulfilled, and, in addition, simple thermodynamic expectations in the limit T → 0 are verified. In fact, in this limit s → 0 according to the third law of thermodynamics and, as ρ B → 0 too, the satisfaction of the condition s/ρ B = const. is insured. We recall (Section V) that, at T = 0, we are in the presence of droplets (states in mechanical equilibrium with the vacuum state (ρ B = 0) at P = 0).
At T = 0, in the first order line, the behavior we find is somewhat different from those claimed by other authors [114,117] where a phenomena of focusing of trajectories towards the CEP is observed. We see that the isentropic lines with s/ρ B = 1, ..., 6 come from the region of symmetry partially restored and attain directly the phase transition; the trajectory s/ρ B = 1 goes along with the phase transition as T decreases until it reaches T = 0; and the other trajectories enter the hadronic phase where the symmetry is still broken and, after that, also converge to the horizontal axes (T = 0). Consequently, even without reheating in the mixed phase as verified in the "zigzag" shape of [113][114][115]118], all isentropic trajectories directly terminate in the end of the first order transition line at T = 0.
In the crossover region the behavior of the isentropic lines is qualitatively similar to the one obtained in lattice calculations [113] or in some models [114,115,119]. The trajectories with s/ρ B > 6 go directly to the crossover region and display a smooth behavior, although those that pass in the neighborhood of the CEP show a slightly kink behavior.
In conclusion, all the trajectories directly terminate in the same point of the horizontal axes at T = 0. As already pointed out in [112], the picture provided here is a natural result in these type of quark models with no change in the number of degrees of freedom of the system in the two phases. As the temperature decreases a first order phase transition occurs, the latent heat increases and the formation of the mixed phase is thermodynamically favored.
We point out again that, in the limit T → 0, it is verified that s → 0 and ρ B → 0, as it should be. This behavior is in contrast to the one reported in [112] (see right panel of Fig.   9 therein), where the NJL model in the SU(2) sector is used. The difference is due to our more convenient choice of the model parameters, mainly a lower value of the cutoff. This can be explained by the presence of droplets at T = 0 whose stability is quite sensitive to the choice of the model parameters. A first point to be noticed is that in the PNJL model, contrarily to what happens in the chiral limit only for the SU(2) sector (m u = m d = 0, m s = 0) where the TCP is found [76], when the total chiral limit is considered (m u = m d = m s = 0), the phase diagram does not exhibit a TCP: chiral symmetry is restored via a first order transition for all baryonic chemical potentials and temperatures (see left panel of Figure 8). Both situations are in agreement with what is expected: the chiral phase transition at the chiral limit is of second order for N f = 2 and first order for N f ≥ 3 [121].
To study the influence of strangeness on the location of the critical points, we vary the current quark mass m s , keeping the SU(2) sector in the chiral limit and the other model parameters fixed. The phase diagram is presented in Figure 8 [122]. The value for m crit s is a subject of debate; those found in lattice [123] or in model calculations [122,124] are lower than the physical strange current quark mass (m s ≈ 150 MeV). We found m crit s ≈ 9 MeV in our model, lower than lattice values [123] and half of the value obtained in NJL model (m crit s = 18.3 MeV [120]), but still consistent with other models of this type [124]. When m s ≥ m crit s , at µ B = 0, the transition is of the second order and, as µ B increases, the line of the second order phase transition will end in a first order line at the TCP. Several TCPs are plotted for different values of m s in the right panel of Figure 8. As already referred, the location and even the existence of the CEP in the phase diagram is a matter of debate [51]. While different lattice calculations predict the existence of a CEP [125], the absence of the CEP in the phase diagram was seen in recent lattice QCD results [79][80][81], where the first order phase transition region near µ B = 0 shrinks in the quark mass and µ B space when µ B is increased [79][80][81]. Due to the importance of the U A (1) anomaly, already emphasized in Section II, and its influence on several observables, it is demanding to investigate possible changes in the location of the CEP in the (T, µ B ) plane when the anomaly strength is modified. In Figure 9 we show the location of the CEP for several values of g D compared to the results for g D 0 , the value used for the vacuum. As already pointed out by K. Fukushima in [70], we also observe that the location of the CEP depends on the value of g D and, when g D is about 50% of its value in the vacuum, the QCD critical point disappears from the phase diagram; the first order region becomes narrower with the decreasing of the g D value. The results show that, in the framework of this model, the existence or not of the CEP is determined by the strength of the anomaly coupling, the CEP getting closer to the µ B axis as g D decreases.

VIII. SUSCEPTIBILITIES AND CRITICAL BEHAVIOR IN THE VICINITY OF THE CEP
In the last years, the phenomenological relevance of fluctuations in the finite temperature and chemical potential around the CEP/TCP of QCD has been attracting the attention of several authors [126]. As a matter of fact, fluctuations are supposed to represent signatures of phase transitions of strongly interacting matter. In particular, the quark number susceptibility plays a role in the calculation of event-by-event fluctuations of conserved quantities such as the net baryon number. Across the quark hadron phase transition they are expected to become large; that can be interpreted as an indication for a critical behavior. We also remember the important role of the second derivative of the pressure for second order points like the CEP.
The grand canonical potential (or the pressure) contains the relevant information on thermodynamic bulk properties of a medium. Susceptibilities, being second order derivatives of the pressure in both chemical potential and temperature directions, are related to those fluctuations. The relevance of these physical observables is related with the size of the critical region around the CEP which can be found by calculating the specific heat, the baryon number susceptibility, and their critical behaviors. The size of this critical region is important for future searches of the CEP in heavy ion-collisions [114].
The way to estimate the critical region around the CEP is to calculate the dimensionless ratio χ B /χ f ree B , where χ f ree B is the chiral susceptibility of a free massless quark gas. Figure 10 shows a contour plot for two fixed ratios χ B /χ f ree B = 1.0 and 2.0 in the phase diagram around the CEP. In the direction parallel to the first order transition line and to the crossover, it can be seen an elongation of the region where χ B is enhanced, indicating that the critical region is heavily stretched in that direction. It means that the divergence of the correlation length at the CEP affects the phase diagram quite far from the CEP and a careful analysis including effects beyond the mean-field needs to be done [127].
One of the main effects of the Polyakov loop is to shorten the temperature range where the crossover occurs [44]. On the other hand, this behavior is boosted by the choice of the regularization (Λ → ∞) [76]. The combination of both effects results in higher baryonic susceptibilities even far from the CEP when compared with the NJL model [128]. This effect of the Polyakov loop is driven by the fact that the one-and two-quark Boltzmann factors are controlled by a factor proportional to Φ: at small temperature Φ ≃ 0 results in a suppression of these contributions (see Equation (19)) leading to a partial restoration of the color symmetry. Indeed, the fact that only the 3−quark Boltzmann factors e 3βEp contribute to the thermodynamical potential at low temperature, may be interpreted as the production of a thermal bath containing only colorless 3-quark contributions. When the temperature increases, Φ goes quickly to 1 resulting in a (partial) restoration of the chiral symmetry occurring in a shorter temperature range. The crossover taking place in a smaller T range can be interpreted as a crossover transition closest to a second order one. This "faster" crossover may explain the elongation of the critical region giving raise to a greater correlation length even far from the CEP. Now we show in Figure 11 the behavior of the baryon number susceptibility, χ B , and the specific heat, C, at the CEP, which is in accordance with [120,126,129].
The baryon number susceptibility is plotted (left panel) as a function of the baryon chemical potential for three values of the temperature. The divergent behavior of χ B at the CEP is an indication of the second order phase transition at this point. The curve for T > T CEP corresponds to the crossover and the other to the first order transition. Now we will pay attention to the specific heat (Equation (27)) which is plotted as a function of the temperature around the CEP (right panel). The divergent behavior is the signal of the location of the CEP. The curve µ B < µ CEP B corresponds to the crossover and the other to the first order transition.
The large enhancement of the baryon number susceptibility and the specific heat at the CEP may be used as a signal of the existence and identification of phase transitions in the quark matter.

IX. CONCLUSIONS
The concept of symmetries is a very important topic in physics. It has given a fruitful insight into the relationships between different areas, and has contributed to the unification of several phenomena, as shown in recent achievements in nuclear and hadronic physics.
Critical phenomena in hot QCD have been studied in the framework of the PNJL model, as an important issue to determine the order of the chiral phase transition as a function of the temperature and the quark chemical potential. Symmetry arguments show that the phase transition should be a first order one in the chiral limit (m u = m d = m s = 0). Working out of the chiral limit, at which both chiral and center symmetries are explicitly broken, a CEP which separates first and crossover lines is found, and the corresponding order parameters are analyzed.
The sets of parameters used is compatible with the formation of stable droplets at zero temperature, insuring the satisfaction of important thermodynamic expectations like the Nernst principle. Other important role is played by the regularization procedure which, by allowing high momentum quark states, is essential to obtain the required increase of extensive thermodynamic quantities, insuring the convergence to the Stefan-Boltzmann (SB) limit of QCD. In this context the gluonic degrees of freedom also play a special role.
We also discussed the effect of the U A (1) axial symmetry both at zero and at finite temperature. We analyzed the effect of the anomalous coupling strength in the location of the CEP. We proved that the location of the CEP depends on the value of g D and, when g D is about 50% of its value in the vacuum, the QCD critical point disappears from the phase diagram. One expects that, above a certain critical temperature T ef f , the chiral and axial symmetries will be effectively restored. The behavior of some given observables signals the effective restoration of these symmetries: for instance, the topological susceptibility vanishes.
The successful comparison with lattice results shows that the model calculation provides a convenient tool to obtain information for systems from zero to non-zero chemical potential which is of particular importance for the knowledge of the equation of state of hot and dense matter. Although the results here presented relies on the chiral/deconfinement phase transition, the relevant physics involved is also useful to understand other phase transitions sharing similar features.