Quark and Polyakov-loop correlations in effective models at zero and non-vanishing density

The aim of this work is to shed light on some lesser known aspects of Polyakov-loop--extended chiral models (namely the Polyakov-loop extended Nambu--Jona-Lasinio and Quark-Meson models), especially on the correlation of the quark sector with the Polyakov loop. We show that the ordering of chiral and Polyakov-loop transitions and their difference in temperature as seen in lattice QCD calculations could be realized with a critical scale of the Polyakov-loop potential that is larger than the one in pure gauge theory. The comparison of the results for the Polyakov-loop susceptibility obtained using the self-consistent medium-dependent quark mass with those obtained while keeping these masses at a fixed value allows to disentangle chiral-symmetry restoration and center-symmetry breaking effects. Furthermore, a confined chirally restored phase is identified by a plateau in the quark contribution to thermodynamics and by sigma and pion spectral functions that coincide but have a small width. We also discuss that, for some large chemical potential values, the explicit center-symmetry breaking is so strong that statistical deconfinement is realized at infinitely small temperatures. Both the missing sensitivity of the Polyakov loop to the quark mass, except at close to the chiral transition, and the Polyakov loop being zero at zero temperature at all chemical potentials, can be interpreted as indications of a missing mechanism which accounts for the quark back-reaction on the Polyakov loop.


I. INTRODUCTION
Recently and in the near future, there is and will be big experimental effort to explore the phase diagram of strongly-interacting matter be it by the energy scan program at RHIC [1], the NA61/SHINE experiment at CERN [2] and or at the future facilities NICA at JINR [3] and FAIR at the GSI site [4]. These laboratory experiments at low and intermediate baryon densities are complemented by the detection of gravitational waves of inspiraling neutron stars that allow to learn about the nature of QCD matter at the largest densities [5]. All these experimental measurements require a theoretical counterpart to interpret and analyse their results. First principle calculations are not yet able to be this match in the medium and large baryon density region. Calculations on a discretised space-time lattice face the infamous sign problem [6] and methods which circumvent it as complex Langevin dynamics are still limited to fundamental theoretical investigations but are not yet connected to phenomenology [7,8]. First-principle continuum calculations using the Functional Renormalization Group, perturbation theory or variational approaches are still in progress towards the true number of quark flavours and quark masses [9][10][11][12].
Therefore, widely used alternatives are frameworks which are based on chiral symmetry and center symmetry and eventually scale symmetry. These symmetries of the QCD Lagrangian are related to two fundamental properties of strongly interacting matter, namely to con-stituent quark masses and confinement in the hadronic phase and the liberation of light quarks in the transition to the quark-gluon plasma. The interaction between constituent quarks that gives them their mass due to spontaneous breaking of chiral symmetry can be described in two ways, either as a point interaction or by exchange of a meson. The former leads to what is called the Nambu-Jona-Lasinio (NJL) model [13] and the latter to the Quark-Meson (QM) model [14][15][16]. Extended by the order parameter of center symmetry breaking, the Polyakov loop, these frameworks allow to explore the phase diagram of strongly-interacting matter [17,18].
Even though the Lagrangian of such models itself is invariant under chiral transformations for massless quarks, the appearance of a non-vanishing quark condensate qq breaks chiral symmetry spontaneously. Therefore, the quark condensate qq is a order parameter for chiral symmetry breaking. The relation between quark masses, chiral symmetry and quark condensate can be exploited to explain the generation of constituent quark masses by spontaneous chiral symmetry breaking. The constituent quark mass of up and down quarks that are confined in protons or neutrons is of the order of one third of the mass of these nucleons, O ∼ 300 MeV which is significantly larger than their current quark masses, O ∼ MeV. While the non-zero current quark masses are responsible for explicit chiral symmetry breaking the constituent quark masses are dynamically generated by spontaneous chiral symmetry breaking, m ∼ qq . With increasing temperature and/or density of the quarks the dynamic chiral-symmetry breaking becomes restored. In the chiral limit of vanishing current quark masses this is a transition of second order which turns at a tricritical point into a first order one. With the explicit chiral-symmetry breaking when the non-zero current quark masses are taken into account the precise second order phase transition becomes washed out into a crossover that can turn in a critical point to a first order transition. The chiral phase structure of these models is for example discussed in Refs. [19][20][21].
The Polyakov-loop field Φ( x) is the appropriate order parameter to study the SU(N c ) phase structure and is associated with Z Nc , center of SU(N c ) [22]. In the pure gauge sector, the phase transition that occurs is related to the deconfinement of colour at high temperature such that the Polyakov loop is an order parameter for the deconfinement transition. Coupling dynamical quarks to the Polyakov-loop field and adding a potential that represents the effective glue potential at finite temperature to drive spontaneous center symmetry breaking gives also information on the deconfinement of quarks. But in this way quarks are not confined in the strict sense since dynamical quarks add an explicit breaking of center symmetry such that the contributions coming from one and two (anti-)quarks are suppressed below the transition temperature but are not vanishing [23].
In this work we will discuss several aspects of the correlation between quarks and the Polyakov loop and investigate how to disentangle effects from restoration of chiral symmetry and from breaking of center symmetry to have a deeper knowledge on how both are correlated in the results of the model. This allows to better understand the current status of comparing results of these kind of models with those of lattice QCD calculations and to shed further light to certain less known results of the model.
In the next section we will give a summarised introduction to the Polyakov-loop-extended Nambu-Jona-Lasinio (PNJL) model, putting special emphasis on the correlation between quarks and the Polyakov loop. In Sec. IV we discuss the correlation of the transition scales of restoration of chiral symmetry and of breaking of center symmetry and how it could be possible to reproduce the splitting that is seen in lattice QCD calculations at vanishing densities. Section V contains the analysis of the effect of the kinetic contribution of quarks to the correlation between quarks and Polyakov loop by comparing results for medium dependent and constant quark mass. The analysis of the effects that can be found in such an investigation are extended in Secs. VI and VII. In Sec. VIII we analyse and discuss the correlation between the quark contribution and Polyakov loop by the chemical potential. Section IX discusses the combination of both effects which is completed by an analysis of a confined chirally restored phase in Sec. X.
We will present results of the model for two flavors. While including the strange quark is important for quantitative predictions, focusing on two light flavours allows us to disentangle non-trivial correlations of the ingredients of the framework giving a more comprehensive view of the physical phenomenon.
II. THE PNJL MODEL IN SU (2) In this section we introduce the Polyakov-loopextended Nambu-Jona-Lasinio (PNJL) model to study the thermodynamics of QCD at finite density and temperature, in the grand canonical ensemble.
The lack of confinement in the NJL models does not allow the study of the very important transition of deconfinement at non-zero temperature. As explained in Sec. I the Polyakov loop can be used as an order parameter for this transition. In this section we will introduce it in the NJL model and describe its features.
Throughout this section we will in particular discuss that if the Polyakov loop brings some sort of confinement ("statistical confinement"), it is not a real confinement. The Fock space structure of this model still contains quark degrees of freedom and only the quark occupation numbers will be modified. As a result when we discuss confinement, we will often not be very rigorous on the vocabulary.
A. Pure gauge sector at finite temperature In the pure gauge sector, the phase transition that occurs is related to the deconfinement of color at high temperature. Following the arguments given in [17,[24][25][26] to study the SU(N c ) phase structure, the appropriate order parameter (associated with Z Nc , center of SU(N c )) is the Polyakov loop field Φ( x): 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 S has been absorbed, P denotes path ordering and the usual notation for the thermal expectation value has been introduced with β = 1/T with the Boltzmann constant set to one (k B ≡ 1). An effective potential that respects the Z 3 symmetry of the original Lagrangian may be built and it is conveniently chosen to reproduce results obtained in lattice calculations. In this approximation, the Polyakov loop field Φ( x) is simply set to be equal to its expectation value Φ =const., which minimizes the potential where T 0 = 270 MeV is the critical temperature for the deconfinement phase transition according to pure gauge lattice findings [27]. A fit of the coefficients a i , b i has been performed in Ref. [28] to reproduce the lattice data [29] for thermodynamics in the pure gauge sector: a 0 = 6.75, a 1 = −1.95, a 2 = 2.625, a 3 = −7.44, b 3 = 0.75, b 4 = 7.5. By minimizing this potential, it was possible to compute the Polyakov-loop expectation value in good agreement with lattice gauge findings [27].
A popular alternative to this potential is the logarithmic one [30] that reads: where The parameters of the effective potential U Log are given by a 0 = 3.51, a 1 = −2.47, a 2 = 15.2 and b 3 = −1.75.
This effective potential with these parameters exhibits the feature of a phase transition from color confinement to color deconfinement through a stronger first-order transition than U P oly with the parameters of Ref. [28], see e.g. Fig. 8 of Ref. [31].

B. Coupling between quarks and the gauge sector
The PNJL model aims at describing in a simple way two of the characteristic phenomena of QCD, namely deconfinement and chiral symmetry breaking. We start from the two flavor NJL description of quarks (global SU c (3) symmetric and chirally invariant point-like interaction), coupled in a minimal way to the Polyakov loop via the following Lagrangian [28] (the range of applicability of this model is limited to temperatures T ≤ 2.5 T c ): where q = (q u , q d ) are the quark fields and where the covariant derivative reads The two-flavor current quark mass matrix ism 0 = diag(m u , m d ) which breaks chiral symmetry explicitly (we work in the isospin symmetric limit and consequently m u = m d = m 0 ). G S is the coupling strength of the chirally symmetric four-fermion interaction. Finally, also a tridimensional ultraviolet cut-off Λ will be introduced in Eq. (9) in order to regularize divergent integrals.
This minimal coupling between quarks and the Polyakov loop is a first way to take into account backreactions of the quark sector to the gluonic sector.
The parameters of the pure NJL sector in Eq. (8) are fixed at zero temperature and density and have the following values: Λ = 651 MeV, G S = 5.04 GeV −2 and m 0 = 5.5 MeV. They are chosen to reproduce the mass of the pion m π = 139.3 MeV and its decay constant f π = 92.3 MeV (obtained in a Hartree + Ring (RPA) calculation) as well as the chiral condensate (order parameter of the chiral symmetry) | qq | 1/3 = 251 MeV. The constituent (Hartree) quark mass is m = 325 MeV. With M N the nucleon mass this means m M N /3. The PNJL model gives a simple explanation of the nucleon as being a state composed of three quarks but their mass is the corresponding dynamical mass. Even with m 0 = 0, m has the same magnitude: essentially the baryonic mass is due to the "glue", the interaction energy carried out by the gauge fields. The mechanism of mass generation via the spontaneous symmetry breaking was one of the first interests of the NJL model. The usual techniques can be used to obtain the PNJL grand canonical potential from the Hartree propagator (see for instance Refs. [23,28,32]), In the above equation E p = p 2 + m 2 is the Hartree single quasi-particle energy (which includes the constituent or dynamical quark mass m and not the current mass The solutions of the mean field equations are obtained by minimizing the grand potential with respect to m, Φ andΦ, namely ∂Ω ∂Φ = 0, ∂Ω ∂Φ = 0 and ∂Ω ∂m = 0. The latter equation coincides with the gap equation, where the modified Fermi-Dirac distribution functions f + Φ and f − Φ have been introduced: The equations presented above, which were introduced for the first time in [23], allow to straightforwardly generalize the results for thermodynamics of the NJL model to those of the PNJL model by replacing the usual Fermi-Dirac occupation numbers by the modified ones given by Eqs. (13) and (14).
The explicit form of the mean field equations for Φ and Φ, which can be found for example in [23], will be given in Sec. II C 2 where we will discuss in detail the correlation between quarks and gauge fields.
Here, we want to point out that we choose not to use the Fock term because, if added, we would obtain the same equation with the replacement G S → G S (1+4/N c ). Since G S is a parameter to be fixed, it is not very important to include the Fock term. We can also notice that the Hartree term is of order O(N 1 c ), the Fock term is O(N 0 c ) and the ring approximation can be shown to be also O(N 0 c ). As discussed in Ref. [33] the Hartree term may be seen as the first term in a 1/N c expansion. Then, for coherence the Fock term should be added when going beyond the mean field.

The influence of quarks on the gauge fields
To perform the study of the influence of quarks on the gauge fields we will push the model in later sections to a very high density scale, probably well beyond the range of validity of the model. However, we find these calculations interesting as they lead us to understand what ingredients are probably missing in the model to get a more faithful representation of QCD in the high density region.
Obviously, in QCD quark fields act upon gauge field via quark loops in the gluon self energy. Such an effect is lacking in the NJL model since there are no gluon fields (although, in some sense gluons are present in the model via the NJL contact interaction as we mentioned earlier).
In the PNJL model there is still no dynamical gluons but there is an interaction between the static gauge field and the quarks via the covariant derivative D µ = ∂ µ − iA µ . As a result, the mean field equation for Φ does have a dependence on quarks (see previous section). Explicitly for two flavors, we have for Φ (the discussion forΦ goes along the same lines): with z + Φ (E p ) and z − Φ (E p ) given by Eqs. (10) and (11) respectively.
The second term in Eq. (15) adds some quark corrections and interesting information can be extracted from this equation as we will discuss later. Without quarks the solution would be the pure gauge one, i.e., (∂U/∂Φ = 0). The model would then have a first-order phase transition.

Statistical confinement
In the modified Fermi-Dirac functions, Eqs. (13) and (14) entering to the grand potential, Eq. (9), the confinement mechanism that exists in Polyakov models can be seen. It is not a true confinement (the quarks are still asymptotically free and we will see that they can be produced e.g. by the sigma meson decay in vacuum) but has been dubbed "statistical confinement".
Indeed, it can be seen that in the grand potential, Eq. (9), the contributions coming from one and two (anti-)quarks are suppressed below T c (when Φ,Φ → 0, confined phase of the model) due to their coupling with Φ andΦ but the 3-(anti-)quarks Boltzmann factor is not. The interpretation is that there are still unconfined (anti-)quarks in the vacuum part (the Fermi sea) of the grand potential but in the thermal bath only 3-(anti-)quarks contributions (a reminder of the fact that in QCD only colorless combinations can exist in the confined phase) are present. This reduces significantly the number of (anti-)quarks in the thermal bath since it requires three times more energy for the (anti-)quarks to be thermodynamically active.
It is known that, at a given value of T and µ, the pure NJL model always overestimates the density (see Ref. [28]), even if it merges for large temperatures with the PNJL model (when Φ → 1). At fixed values of T and µ, the PNJL value for n q is much lower than in the NJL case. In fact, all the possible contributions to the latter turn out to be somehow suppressed: the one-and two-quark contributions because of Φ,Φ → 0, while the thermal excitation of three quarks has a negligible Boltzmann factor. We would be tempted to identify these clusters of three dressed (anti-)quarks with precursors of (anti-)baryons but no binding for these structures is provided by the model. In any case it is encouraging that coupling the NJL Lagrangian (which parameters are chosen to reproduce zero temperature properties) with the Polyakov-loop field (described by a pure gauge effective potential) leads to results pointing into the right direction.
Another important impact of statistical confinement is that the transition will occur in a shorter range of temperature than in the NJL model. It is seen from Eqs. (10) and (11) that before the transition, the pressure (for example) is kept low by the effect just discussed. Then, quark degrees of freedom are liberated in a short temperature range when Φ → 1. As seen in Ref. [34] the transition in the PNJL model is indeed much "faster" that the one in the pure NJL model.
The mesonic contribution to the pressure is also important (about 10%) and this effect has to be included in the mean-field description (e.g. with the Beth-Uhlenbeck formalism [35] which is another way of doing RPA) or to a model with mesonic degrees of freedom such as the Polyakov-loop-extended Quark-Meson (PQM) model [36]. This allows to fit nowadays well the most recent lattice results with physical quark masses. A detailed discussion of the results in mean-field approximation can be found in [34].

III. COINCIDENCE BETWEEN CHIRAL AND DECONFINEMENT TRANSITION
One of the first successes of the PNJL model was the observation that without any additional tuning, chiral and deconfinement transitions almost coincide in a small range of temperatures at zero density (see Fig. 7). It was an interesting test as LQCD calculations confirmed this coincidence (even if recent lattice results show there is no perfect coincidence [37]).
What is not obvious is that by mixing the deconfinement scale T 0 = 270 MeV and the chiral restoration scale (of about 220 MeV in the NJL model) the two transitions automatically coincide at a lower temperature. The PNJL minimal coupling is enough to have some sort of back-reaction effect that produces this matching. This is also a quite stable feature [28]: we can take T 0 = 190 MeV and still get a quite good coincidence. The motivation to change the value of T 0 is to get a better agreement with the value of the LQCD chiral transition temperature. As pointed out, it can be justified since in an effective model we are mixing physical sec-tors with different scales in a symmetry based (Landau) framework. It is then understandable that the absolute scale of the two sectors can be slightly adjusted to get a better agreement with the phenomenology. In other words, we authorize ourselves to consider T 0 as a free parameter of the model with a loose constraint on it coming from LQCD, exactly as it can be done in the chiral sector with the constituent quark masses despite that there is some estimate of these from measurements.

IV. ADJUSTING CHIRAL AND DECONFINEMENT PHASE SCALES
The most problematic aspect for PNJL type models is probably the fact that the difference between the chiral transition and the raise of Φ in the QGP in recent lattice data is larger than previously seen, see Fig. 1. To show this, we use a recent PQM model [36] because we need a good agreement between the model and LQCD data. In this figure, t is the reduced temperature adjusted to the chiral crossover temperature T c : t = T /T c − 1. The pressure plot, or the quark number density one, shows that at T c (t = 0) quark degrees of freedom are liberated. But there are two mechanisms responsible for this liberation: the chiral restoration and the deconfinement transition. This can be easily studied in effective models since we have a control on how these phenomena are coupled.
In Fig. 1, bottom left, it is possible to see that at t = 0 chiral symmetry goes toward its restoration. There is a liberation of thermodynamical degrees of freedom simply because m → m 0 hence the quark Boltzmann factor has a bigger contribution to the pressure. This occurs at the chiral transition scale. On the bottom right panel, it is seen that the results for Φ from LQCD and PQM model calculations do not agree. The (statistical) deconfinement occurs at a higher temperature than in the PQM model (at t = 0, the PQM model gives Φ 0.5 whereas it is only 0.1 for LQCD results). I.e. the (statistical) deconfinement scales differ in LQCD data and in the PQM model.
One subtlety concerning the Polyakov loop is that there are different order parameters. The standard order parameter is the expectation value Φ [A 0 ] . The functional dependence indicates that the Polyakov loop derives from the temporal component of the gauge field A µ . The expectation value of the later A 0 relates to another Polyakov-loop order parameter, The LQCD results are for Φ since this is the easily accessible quantity in these results. But in Polyakov-loop extended effective models for QCD the Polyakov-loop potential is adjusted to lattice results on Φ in pure gauge theory while it is actually the gauge field that appears in the fermionic determinant initially, see Eq. (9) which is then rewritten into a dependence on the Polyakov loop Φ in Eqs. (10) and (11). Both Polyakov-loop order parameters are related by Jensen's inequality in the following way: [40,41]. This condition shows that the transition scales of both Polyakov-loop order parameters might differ with that of Φ being larger than that of Φ [ A 0 ]. The results shown in Fig. 1 are consistent considering this condition.
The origin of these both Polyakov-loop order parameters, their derivations and their relation is discussed in detail in Ref. [42]. The temperature dependence of Φ [ A 0 ] is calculated in continuum approaches such as the Functional Renormalization Group, with Dyson-Schwinger equations and the 2PI formalism as well as in the Hamiltonian approach and covariant variational approach and in perturbative approaches. For pure gauge theory, available results on Φ [ A 0 ] in different continuum approaches are shown together with those for Φ [A 0 ] of different lattice calculations in Fig. 2. Interestingly enough, for the perturbation theory calculations, the leading order result agrees with the results of the other continuum approaches while evaluating the next-to-leading-order result with the leading-order relation between expectation value of the gauge field and that of the Polyakov loop (see also Ref. [47]), it is very close to the lattice results.
The information given by LQCD calculations is that the scales of chiral symmetry restoration and center sym-  metry breaking do not match perfectly. In general, in the PNJL/PQM model families, this scale separation is not achieved or worse, with common parametrization as the one we use in PNJL (T 0 = 190 MeV), the scale hierarchy is inverted, see Fig. 7.
In effective models, one can have control over this hierarchy as we show in the following. The shortcoming however is that the output of the model for the absolute temperature of the phase transition will not be correct. This is a problem but not to the point that the model becomes meaningless. Indeed, in the spirit of Landau theory, two systems with the same universal behavior cannot be directly compared but they have to be compared relatively to a given scale. This is of course not a complete justification for what we will do in the following but it will illustrate which information is possible to obtain about some mechanism that is probably missing in the model. We recall that a scale adjustment was already done at the beginning of the PNJL model: for example, in Ref. [28] the T 0 scale is modified to get a better absolute agreement with the chiral transition temperature given by LQCD results.
So the point is that we want to have in the model the correct scale hierarchy for the transitions. In the PNJL model we have two sectors to determine the respective scales; the NJL one (fitted to chiral symmetry breaking phenomenology in vacuum, a scale related to the strength of the condensate) and the gluonic one (fixed by pure gauge results at finite temperature with the scale T 0 ). The coupling is done via the covariant derivative but if we allow to vary T 0 then we can control the relative scales of the transitions. This is done in Fig. 3: we see that the correct hierarchy can indeed be achieved with T 0 400 MeV but the price to pay is that both chiral and deconfinement transitions occur at too high temperatures.
To get a better phenomenology one could compute the PNJL properties not relatively to an absolute scale but relatively to t where T c is the chiral transition temperature. Then, one adjusts the NJL parameters to get a better mesonic phenomenology when expressed relatively to this scale. The latter adjustment is needed because without it, the pion mass would be too low. Indeed, if we have to rescale everything to get the correct absolute scale, it would be needed a correction of about a half to send T c 400 MeV to the correct T c 200 MeV hence the pion mass would be around 65 MeV.
We are not sure if such a full reparametrization is possible but anyway the message is that some mechanism is missing to get the correct scale hierarchy and we should try to understand which mechanism would allow to have the same effect than our rather artificial increase of T 0 . Since up to now no model can reproduce this fact, we probably need to look outside of the existing mechanisms of the PNJL/PQM models. Our guess is that introducing a minimal dynamics to the gauge sector [48] may be the correct way.
We also notice that in the PQM model calculation shown in Fig. 1 the pressure and quark density are probably a little bit overestimated for t > 0 when compared to LQCD. If the correct scale hierarchy could be imple- mented, maybe the later liberation of quarks via statistical deconfinement could reduce this excess. This is an a bit far fetched conclusion since the excess is low (statistically it seems around a small but constant 1 sigma deviation) and in PQM model the mesonic degrees of freedom tend to survive in the QGP phase [49] Finally, for what concerns this aspect of quark and Polyakov loop correlations, we also notice in Fig. 3 that the full calculation for Φ from Eq. (15) and the one keeping the mass m at a fixed value in Eq. (15) coincide, except when the difference between the chiral temperature and the fixed m deconfinement temperature is smaller than about 25 MeV. This means that the order parameters Φ and the chiral condensate are weakly coupled by the mean field equations because except for a small 25 MeV Λ QCD scale difference, the decorrelated calculation (fixed m) gives the same result as the correlated one. We will discuss this kind of analysis in more detail in the next section.
For the quarks, at low T 0 the chiral transition is below its NJL value (around 220 MeV) and it happens at a large (but below 1) value of the Polyakov loop while at high T 0 the transition occurs for small value (around zero) of the Polyakov loop. Thus, in the PNJL model the chiral transition always behave differently from what happens in the NJL model due to the action of the Polyakov loop, explaining the success of the model. The shortcoming of the model is essentially in the behavior of the Polyakov with respect to quarks but the converse seems to be fine.
The general conclusion is that quite certainly this class of models underestimate the correlation for the lack of dynamical quark loop effect on the Polyakov loop. The only effect present is due to the presence of thermodynamical quarks (via the quark Boltzmann factor in the mean field equation) that are due to the minimal couplinḡ qA µ q. But since A µ is not dynamical this term will not generate quark loops in a diagrammatic approach (A µ acts as an external constant gauge field only).

V. DYNAMICAL THERMAL QUARKS EFFECT
In this Section we will consider the effect of the kinetic term of quarks in Eq. (15). At zero density (µ = 0), and with infinite quark mass, the second term disappears as it should. Indeed, quarks are then no more dynamical (they are quenched) and the system is in its pure gauge or Yang-Mills limit. The Z 3 symmetry is then exact (for the grand potential) but spontaneously broken at T > T 0 .
When quarks are considered dynamical again, by lowering their mass approaching the respective physical values, the Z 3 symmetry gets also an explicit breaking term Sec. 2.2.2 and 2.3.2 in Ref. [36]) because of the kinetic term in the Lagrangian. This mechanism can be understood in detail by solving the mean field equations, not self-consistently as usual, but simply by fixing in Eq. (15) the mass m to a given value. In Figs. 4 and 5 we plot the behavior of the Polyakov loop for different masses, and hence changing the ratio between the mass energy E m and the kinetic energy E k of a quark: increasing the mass reduces the ratio E k /E m ; when this ratio is zero (quenched quarks) the symmetry is restored. To better illustrate if there is a stronger or weaker explicit breaking of the Z 3 symmetry, we use as Polyakov-loop potential a parametrization that has a strong first-order transition in the pure gauge sector, (the U Log potential, Eq. (5)). If the transition becomes smoother (crossover like transition) it indicates that we have a stronger explicit breaking of Z 3 symmetry.
These figures also show the impact of varying the chemical potential. Indeed, there is a contribution to the kinetic energy due to the Fermi momentum of quarks: increasing µ increases the ratio E k /E m and the symmetry is more broken. We will discuss further this effect in Sec. VII but for now, we focus on the zero density situation where there is only one contribution to the kinetic energy.
We observe that a lowering of the mass term introduces a stronger explicit breaking of Z 3 symmetry as seen by the large crossover obtained and the respective smooth Polyakov loop susceptibility. We also notice that for m below the scale Λ QCD the qualitative behavior of the Polyakov loop becomes relatively insensitive to the value of m. In Fig. 6 it is seen better that Φ can be computed with m = m 0 or m = m vac 325 MeV without significant difference.
We notice that, except very close to the chiral transition (where m changes quickly) χ Φ is not much altered by keeping the mass constant. So, even the transition properties are almost not affected by the precise value of the mass. It is probable that this lack of sensitivity to the mass of the quark is a sign of a missing mechanism to take into account the quark back-reaction on Φ because of the lack of dynamical quark loops in its calculation. For this reason, in the literature it is proposed to add phenomenologically more sources of back-reaction by using for example a scalar coupling G that depends on Φ, the entangled PNJL (EPNJL) model [50]. Fitting this dependence is quite difficult (due to the lack of enough data) and some authors proposed to use the Roberge-Weiss transition at imaginary chemical potential [51] to fit this dependence. Other sources of back-reaction may be added by using a µ-dependent gauge transition temperature T 0 (µ) [52] or even a µ-dependent Polyakov potential U(Φ; T, µ) [53].
We will not discuss further the EPNJL model but here we want to stress that, before going to this kind of model with reinforced back-reaction, it is interesting to understand the multiple aspects of back reaction already in the PNJL model. Furthermore, models with many sources of back-reaction (G(Φ), T 0 (µ)) introduce many new parameters and it is not easy to consistently parametrize them without overestimating the back-reaction. In this regard, LQCD data at finite density would be very helpful or high density constraints (e.g. the measurement of the position of the CEP in experiments). It is interesting to note that LQCD data at finite imaginary chemical potential exists [54] and it is believed that the continuation iµ → µ is analytical (it can be seen easily in the PNJL model grand potential) so constraints at imaginary µ are relevant for models to be used at real µ.
It is also worth mentioning a last remark concerning Fig. 6: as already discussed in Sec. II C 3, the main effect of the Polyakov loop onto the quark condensate is to obtain a steeper crossover transition and in Fig. 6 (lower left panel) we see indeed that when χ Φ is at its maximum it also introduces a dominant peak in χ σ .

VI. MEAN FIELD PHASE DIAGRAM
Using the techniques described in [34] we investigate the phase structure of a common two flavours PNJL model, computed in this case with the logarithmic Polyakov-loop potential 1 . To investigate the deconfinement transition we study the Polyakov loop susceptibility, and extract from it the characteristic transition temperature of center symmetry breaking as a sign of deconfinement. Some of the details of the numerical calculations for a similar model, that are not of importance for our quantitative discussions here can be found in [57]. To obtain the susceptibilities, the "easy way" is to compute them as numerical derivatives. This is both inefficient (it is quite long since for each step in temperature the mean field equations have to be solved several times) and inaccurate. The inaccuracy will show up as a numerical fluctuation in the transition line plot which shows the maxima of the susceptibilities. A more accurate and faster way is to compute analytically (or at least semi-analytically in the sense that the non analytic thermal integrals are kept as numerical integrals) all first and second order partial derivatives of the grand potential with respect to the order parameters and T and µ and then combine them to get the desired susceptibilities (see [57,58]).
Results for the order parameters m and Φ are pre-sented in Fig. 7. Concerning the chiral properties, they behave as discussed in the introduction and we observe a first-order transition line that starts at zero temperature and ends with increasing temperature at the critical end point (CEP), where a second-order phase transition occurs, continued by a crossover for further increasing temperatures.
For what concerns the deconfinement transition we see that, except for chemical potentials 250 MeV where χ Φ has several maxima, the order parameter seems to always follow a crossover transition. These multiple maxima are simply due to the presence of m in the mean field equation for Φ: when dm/dT has multiple maxima also has Φ. If we look at the Polyakov loop susceptibility in Fig. 8, we see that there is a peak at the chiral transition temperature but physically this is not the signal of the deconfinement transition. We see that in the region between the different maxima, Φ increases only about 5% and stays below 0.2: this is not a (first-order) deconfinement transition because quarks are still statistically confined when Φ = 0.2. There is also a bump in χ Φ , and there Φ changes between 0.4 and 0.7 and the system reaches the statistically deconfined phase. Hence, we will consider this temperature as the one of the deconfinement crossover.
For sake of completeness, we have plotted in the phase diagram in Fig. 7 besides the maximum of χ Φ and its inflection points (that estimates the "size" of the crossover region), also the minimum of χ Φ , when it exists, between the deconfinement crossover bump and the first-order chi- ral symmetry phase transition peak. It is worth noticing that both transitions seem to be strongly correlated below µ 250 MeV and then they become decorrelated at higher chemical potentials, a sign of the mutual influence between chiral and deconfinement sectors. The topic of the quark -Polyakov loop correlation deserves therefore special attention as discussed in Sec. II C 2.

VII. QUARK CONDENSATE AND POLYAKOV LOOP CORRELATION
We continue our study focusing on finite density phases.
At zero density, the pseudocritical temperatures for partial chiral symmetry restoration and deconfinement are very close, differing by approximately 20 MeV [37].  At finite density the status is not clear (different model parametrizations give very different results). Experimentally, there is a lack of high temperature and density experiments. Furthermore, it is also difficult to extract information about the chiral sector only or the deconfinement sector only. In our model (this feature is model dependent and will eventually disappear by considering a T 0 (µ) parametrization, see Ref. [59]) we will study the high density phases and the possible opening of a new phase near the CEP. Close to the CEP the chiral and deconfinement transitions decouple as seen in Fig. 7.
In Fig. 9 the Polyakov loop susceptibility is computed with the full mean field equations, and with a mass fixed to the mean value of its value in the vacuum and the current quark mass m 0 , denoted m mean (dashed lines). We can see that for a chemical potential that gives a firstorder phase transition for the chiral transition (red full line), there is a peak in the Polyakov loop susceptibility at the chiral transition. Having in mind the discussion in the previous section, at this peak, the Polyakov only varies from 0.2 to 0.3 There is also a bump in dΦ/dT at higher temperature and around this bump, Φ goes from 0.4 (essentially confined phase) to 0.7 (essentially deconfined). Our interpretation is the following: it is not possible to say that a first-order phase transition occurs for Φ even if there is a peak in χ Φ because the value of the Polyakov loop remains small. However, there is a smooth crossover towards the deconfined phase at higher temperatures. These two transitions must be considered and this interpretation opens the possibility to describe an intermediate phase, sometimes called quarkyonic phase, that we prefer to call confined chirally restored phase (CCS).
Indeed, it can be explicitly shown that in this phase quark matter is still (statistically) confined. However, we postpone the detailed discussion to Sec. X.
In Fig. 9 it can be seen that when using a constant mass, the chiral peak in the Polyakov loop susceptibility disappears (dashed lines). This is understandable since we argued that there is not a significant variation in Φ at these temperatures but only a brutal mass change that creates the large derivative. Φ itself does not vary much since we saw that it is not very sensitive to mass variations below Λ QCD .
When the mass is fixed at a chemical potential lower that the one of the CEP, the maximum is approximately the same. This shows that taking m fixed really allows to decorrelate the chiral condensate effects because the shape of the Polyakov loop susceptibility does not change significantly when the mass varies strongly (in this interval the mass goes from 320 MeV to 150 MeV). All the remaining effects are only due to the chemical potential, an aspect that we will discuss in the next section about the Fermi momentum effect.
Another important aspect is that the use of a fixed mass is also very useful numerically. If we want to draw a phase diagram but the exact behavior of the Polyakov loop near the CEP is not so relevant, we can use this trick and get a faster algorithm since only two mean field equations have to be solved 2 . This method can be helpful even to compute the full solution: it can greatly increase the calculation speed if we start to compute the approximated fixed mass value and then incorporate it in the solving algorithm for the full equations. This kind of algorithm (root polishing when the approximate solution is known) is particularly helpful when the dimension of the problem increases (here three dimensions since there are three mean field equations; including the strange sector and without isospin symmetry it would be five dimensions).
It is possible to qualitatively compare both methods, fixed mass and full calculation, by looking at the phase diagram in Fig. 10 (here the Polyakov potential is the Polynomial one but qualitatively our results are independent on the choice of potential ). The conclusion is that except in the region of the CEP it does not differ from the previous one found in Fig. 7, Sec. VI. The extra features will be discussed in the next section since they concern large density effects.

VIII. FERMI MOMENTUM EFFECT
In the previous section, we have seen in detail one important aspect of the chiral transition/deconfinement transition correlation by carefully considering the effect of the quark masses on the field equations. As mentioned, there is another quark effect, in the mean field equation via the quark chemical potential, that originates from a non zero Fermi momentum of the quarks. Because of the covariant derivative and the inclusion of µ in the Lagrangian, there is a dependence on µ in Eq. (15) that can have a strong effect on the system. In fact, the Fermi energy will act exactly as any kinetic energy (in the Lagrangian µγ 0 acts as ∂ 0 γ 0 ) and increases the explicit breaking of Z 3 symmetry. In Fig. 9 we computed χ Φ for several chemical potentials. It is verified that, at some point, the breaking of the symmetry is so strong that the transition cannot be defined even as a crossover anymore (at least according to our rigorous definition ). From the moment on where the Polyakov loop susceptibility does not have a maximum anymore, our interpretation is that the transition goes so fast from Φ = 0 to 1 (and without the characteristic of a crossover transition where correlation lengths quickly increase then decrease) that we can almost say that the system is deconfined even at zero temperature (in Fig. 10, Φ is already 0.5 at T = 50 MeV for large values of µ). Φ is forced to be zero at zero temperature but in the large µ limit, Φ = 1 at any non-vanishing temperature. Hence, the model basically tell us that, apart a small missing ingredient to allow Φ to be non zero at zero temperature, a statistical deconfinement induced only by the density is possible. This missing ingredient is again probably related to quark back-reaction on the gauge fields.
Here we want to point out that it is quite remarkable that, in its simple version, the PNJL model already allows a deconfinement transition caused only by the density (with a small temperature), without any sort of ad hoc back-reaction added.
In the phase diagram presented in Fig. 10 more features can be read. It is interesting to have an idea of the value of the order parameter, Φ, with respect to the transition lines, so we add the Φ = 0.5 characteristic line (pink dashed curve). More importantly, we show in this diagram that there is a kind of "crossover end point" in this model calculation. At very high µ (see also Fig. 9) the explicit breaking of Z 3 symmetry is so strong that χ Φ has no longer a maximum (except at T = 0). Somehow this explicit symmetry breaking is dominant in the Lagrangian and after this point it is no longer possible to talk about a transition, even a crossover one. The symmetry is completely destroyed by this breaking term.
As a conclusion, we emphasize that only the careful assessment of the sources of correlations in the mean field equations allowed us to see those effects. Again, we advocate that even in this simple model some mechanism as the Fermi momentum action can be seen. The use of more realistic models (e.g. SU(3) PNJL model) would also be interesting to get more faithful results. However, we think that those models, with many interplays between physical sectors whose contributions are controlled by parameters adjusted to phenomenology, are complementary to the simpler models and both have their merit.
For example, if we want to answer a qualitative question about a physical phenomenon, like if the existence of a CCS phase will favor or disfavor the existence of compact star with quark matter cores, we believe that the results given by a SU(2) PNJL model are much more useful as a guide, while a quantitative answer must be given by a more realistic model where all physical mechanism contribute.

IX. COMBINED EFFECTS OF CHIRAL CONDENSATE AND FERMI MOMENTUM
To complement our discussion on the influence of changing the chiral condensate (by varying m) and the Fermi momentum (by varying µ) on Φ, Fig. 11 presents the results for the deconfinement temperature obtained when varying both µ and the fixed mass m f ix .
Taking µ = 0, the deconfinement transition for m f ix = 0 is a crossover with a treansition temperature about the same as the chiral one (for this model parametrization it is about 220 MeV). This corresponds to the usual situation of the PNJL model when both transitions almost coincide at zero density. As already explained when discussing Fig. 4, by increasing the mass (thus decreasing the explicit Z 3 symmetry breaking) the temperature increases until a value where the crossover becomes a firstorder phase transition, at a temperature of about 270 MeV (which is the value of T 0 in this parametrization).
Now, looking at m f ix = 0 and increasing µ, we see a lowering of the deconfinement crossover temperature until a value where the crossover cannot be defined anymore. The symmetry breaking induced by the Fermi momentum is so strong that after this point, Z 3 is no more a valid symmetry of the model.
Finally, by varying both µ and m f ix it is always possible to go to a first-order phase transition if the mass is sufficiently high compared to µ in order to counteract the explicit breaking of Z 3 symmetry.

X. CONFINED CHIRALLY RESTORED (CCS) PHASE
In this section we will take a closer look to the CCS phase, a phase where chiral symmetry and deconfinement are decorrelated.
In the PNJL model it is always necessary to remember that quarks are not confined in the strict sense and when we talk about confinement we mean statistical confinement in the sense given by the definition in Sec. II C 3. These type of models cannot tell us the true nature of quark confinement.
To explore the CCS phase we want to count the number of thermal quark degrees of freedom (the only one affected by statistical confinement). For this reason we need to subtract the quarks in the vacuum, present due to the lack of confinement (we have seen that their contribution is quite reasonable by looking to the quark number density in Sec. II C 3). Hence, we calculate a relative pressure difference by subtracting the pressure at zero temperature and normalizing it to the Stefan-Boltzmann pressure: ∆P ≡ P (T, µ) − P (T = 0, µ) P SB (T, µ) − P SB (T = 0, µ) . Relative pressure difference (color map), in the (T, µ) plane, super imposed on the phase diagram of Fig. 7 .
This quantity is adequately related to the thermal quark degrees of freedom and indeed it shows signs of statistical confinement in the CCS phase (see Fig. 12). In the hadronic phase, ∆P is zero (the black area) meaning that the PNJL model will not thermodynamically create any quarks in the hadronic phase (besides those already present in vacuum). In the QGP phase (the yellow region), thermal quarks are very active thermodynamically. Finally, in the CCS phase a quite stable "plateau" at around 0.2 is found. At a given (small) temperature, there is a discontinuous jump from 0 to 0.2 when increasing µ across the first-order transition. So, in the CCS phase some quarks have been liberated by the system due to the chiral transition but their contribution to thermodynamics remains low. It seems that in this phase, where the chiral symmetry is restored, the model tells us that quarks still are far from being asymptotically free. Phenomenologically, this could have rather important effects on experimental observables in HICs (viscosity, hadronic abundances, etc). If this hypothetical phase exists in QCD, the PNJL model cannot say if quarks are really confined (no asymptotic states) or if they are deconfined but the interaction is still strong. Performing the same calculation in the pure NJL model, i.e. without statistical confinement, one observes at small temperature that ∆P has the same jump at the chiral transition and then increases slowly without a "plateau" behavior.
It is also relevant to notice that, at those densities, an important missing effect is the presence of diquark condensates (the strange sector is also important). This point deserves further investigations.
To confirm this finding and that it may lead to significant differences in HICs, we study a more direct observable, namely we check if mesonic probes act as it is expected in a confined but chirally symmetric phase.
We compute the spectral function (see [23]) in the meson to quarks channel at zero momentum and report FIG. 13. Upper panels: Pion and sigma spectral function in the hadronic phase (left) and the QGP phase (right). Confining regime: σ is a narrow resonance; the chiral symmetry is broken and qq = 0; Φ 0. Deconfinement: σ and π spectral functions are almost the same as an uncorrelatedqq pair; the chiral symmetry is restored and qq 0; Φ = 0. Lower panel: Pion and sigma spectral function in the CCS phase. It looks like a non asymptotically free regime (statistically confined at least) since σ and π are still narrow resonances. But chiral symmetry is restored: qq 0; Φ 0. them on Fig. 13. The finite width of the sigma in vacuum indicates that it is not confined (this channel is open) but there is a strong interaction and thus it is a sharp resonance. The expected behavior for mesons in a CCS is observed in Fig.13: chiral symmetry is effectively restored when the sigma and pion spectral functions essentially coincide [23], the statistical confinement is still strong as it can be seen by the fact that their width is small (right upper panel) compared to the full deconfined phase (lower panel).
Finally, we remark that other parametrizations, for example with a chemical potential dependence of T 0 , will essentially close this CCS phase [59]. Yet, we think it is still interesting to try to understand better the CCS phase. We also think that T 0 (µ) is a too strong a priori generalization to introduce in the model. This de-pendence is justified theoretically in Ref. [52], where it is argued that this minimal necessary generalization is needed because, without a µ-dependent b 2 in Eq. (4), the confinement-deconfinement phase transition has a higher critical temperature than the chiral phase transition, being an unphysical scenario since QCD with dynamical massless quarks in the chirally restored phase cannot be confining because the string breaking scale would be zero. However, the problem is that it is badly constrained. By using this parametrization the hypothesis we have is basically that the deconfinement and chiral restoration will coincide Of course the CCS phase will then not exist anymore but in our opinion this is not an argument against CCS existence. It can eventfully be an argument against the CCS phase only if this a priori generalization is correct but which has to be put on more firm basis (in particular better constrained).
In the meantime we think it is still interesting to consider models with a CCS phase even if in the future it will be contradicted by experimental data, the ultimate judge. Anyway, it is still a nice exercise to see what physical information can be gathered by using effective models.

XI. CONCLUSIONS AND OUTLOOK
We started this working by showing the current status on how far effective models based on the restoration of chiral symmetry and the breaking of center symmetry (to describe the deconfinement transition of QCD) are able to reproduce the result of lattice QCD (LQCD) calculations. Based on the observation that there is a substantial quantitative difference in the temperature dependence of the Polyakov loop, regarding the chiral condensate one, we have analysed several aspects of the correlation between the chiral sector and the Polyakov loop. Examining the pseudo-critical temperatures of chiral and Polyakov-loop transition in dependence of the critical scale of the Polyakov-loop potential T 0 we see that the order of chiral and Polyakov-loop transition and their difference in temperature as seen in LQCD calculations would be realised with a T 0 larger than the one in pure gauge theory. The absolute temperature scale of the QCD transition would then be larger than in LQCD but in the spirit of Landau theory all scales could be rescaled to arrive at the correct value for T c .
For studying the effect of the kinetic contribution of quarks to the correlation between quarks and Polyakov loop we compared the behaviour of the later between calculations using a self-consistent medium-dependent quark mass and for several values of constant quark masses. We observe that for values of the dressed quark mass m below the scale Λ QCD the qualitative behaviour of the Polyakov loop becomes relatively insensitive to the exact value of m. Except very close to the chiral transition (where m changes quickly) the Polyakov-loop susceptibility is not much altered by keeping the mass constant so that even the transition properties are almost not affected by the precise value of the mass. We interpret this lack of sensitivity to the mass of the quark as sign of a missing mechanism to take into account the quark back-reaction on the Polyakov loop because of the lack of dynamical quark loops in its calculation. Comparing the Polyakov-loop susceptibility for medium-dependent and constant mass at finite density allows to disentangle a smooth crossover proper to breaking of center symmetry from a peak in the susceptibility introduced purely by a chiral first order transition but at which the value of the Polyakov loop remains small, so in the confined regime. This picture opens the possibility to discuss an intermediate phase, general known as quarkyonic phase, that is in our model a confined chirally restored phase. Probing the number of thermal quark degrees of freedom within this phase, indeed one sees signs of statistical confinement with a quite stable "plateau" where some quarks have been liberated due to the chiral transition but where their contribution to thermodynamics remains low. Furthermore, meson spectral functions behave as it is expected in a confined but chirally symmetric phase: the sigma and pion spectral functions essentially coincide but their width is small. Another contribution to the correlation between quarks and Polyakov loop at non-vanishing density comes from the fact that the Fermi energy acts exactly as any kinetic energy and increases the explicit breaking of center symmetry. From some large chemical potential on this explicit breaking of the symmetry is so strong that the Polyakov-loop susceptibility does not have a maximum anymore except at T = 0. Hence, apart a missing ingredient to allow the Polyakov loop to be non-zero at zero temperature, a statistical deconfinement induced only by the density seems to be realisable within the model. This missing ingredient is again probably related to quark back-reaction on the gauge fields.
Overall, the goal of this work is to shed light to some less known aspects of Polyakov-loop extended chiral models to motivate certain paths for future improvements and detailed quantitative analysis.