Thermodynamical properties of strongly interacting matter in a model with explicit chiral symmetry breaking interactions

We analyse the effects of the light and strange current quark masses on the phase diagram of QCD at finite temperature and vanishing baryonic chemical potential, computing the speed of sound, the trace anomaly of the energy momentum tensor and the fluctuations and correlations of the conserved charges associated to baryonic, electric and strangeness numbers. The framework is a known extension of the three flavor Nambu Jona Lasinio model, which includes the full set of explicit chiral symmetry breaking interactions (ESB) up to the same order in large $N_c$ counting as the 't Hooft flavor mixing terms and eight quark interactions. It is shown that the ESB terms are relevant for the description of a soft region in the system's speed of sound and overall slope behavior of the observables computed. At the same time the role of the 8q interactions gets highlighted. The model extension with the Polyakov loop is considered and the results are compared to lattice QCD data.

We analyse the effects of the light and strange current quark masses on the phase diagram of QCD at finite temperature and vanishing baryonic chemical potential, computing the speed of sound, the trace anomaly of the energy momentum tensor and the fluctuations and correlations of the conserved charges associated to baryonic, electric and strangeness numbers. The framework is a known extension of the three flavor Nambu Jona Lasinio model, which includes the full set of explicit chiral symmetry breaking interactions (ESB) up to the same order in large Nc counting as the 't Hooft flavor mixing terms and eight quark interactions. It is shown that the ESB terms are relevant for the description of a soft region in the system's speed of sound and overall slope behavior of the observables computed. At the same time the role of the 8q interactions gets highlighted. The model extension with the Polyakov loop is considered and the results are compared to lattice QCD data.

I. INTRODUCTION
The study of the thermodynamical properties of strongly interacting matter is one of the open present day theoretical and experimental challenges. On the theoretical side the well known problems of the ab initio approach of lattice QCD (lQCD) when dealing with finite chemical potential as well as the will to grasp a deeper understanding of the interplay of the various underlying mechanisms in strongly interacting matter can be seen as an incentive to the use of moderately complex effective lagrangians. On the other hand at vanishing chemical potential the growing confidence in lQCD results means that the agreement with these is increasingly used as a way to establish the success of other theoretical predictions.
The values of the current quark masses constitute undoubtedly one of the most relevant inputs in the study of the QCD phase diagram, as the quark condensates become not exact order parameters for the chiral phase transition. Model estimates of the size of the condensates at the critical transition points show that they may be significantly larger than the bare ones [1][2][3][4][5][6][7], indicating that non-perturbative effects are still effective in spite of the transition. After the transition, a more or less slow convergence to the bare values of the condensates depends naturally on the size of the current quark masses and how fast the perturbative regime of QCD is reached, where bulk thermodynamic observables are conditioned by the Stefan-Boltzmann limit pertinent to an ideal quark-gluon gas.
The chiral critical end point (CEP) which according to numerous model calculations is expected to occur, separating a region of first order transitions at higher baryon * jmoreira@uc.pt † aaosipov@jinr.ru chemical potential µ B and lower temperatures T from the crossover behavior at lower µ B and higher T, is not yet established. A second order transition is likely to occur at this point in the chiral limit of light quarks and an infinitely heavy strange quark [8]. Model results, using quark masses close to physical values, differ drastically regarding its possible location. The hope that its eventual location may be narrowed down in lQCD, using the reweighting technique to extend lattice calculations from µ B = 0 to finite values [9], should the nature of the transition be second order, requires still a detailed analysis of the quark mass and volume dependence [10], on which hinges the accuracy of Lee-Yang zeros of the lQCD partition function [11].
Meanwhile extensive lQCD studies by different groups report that the algorithmic difficulties that prevented the use of the light physical quark masses have been mostly overcome, as well as spurious taste breaking effects for staggered discretization schemes, making it possible to achieve a realistic hadron spectrum [12,13]. Their findings converge to the by now commonly accepted understanding that along the µ B = 0 line no genuine phase transition occurs. A crossover takes place around T ∼ 155 MeV in a T interval of roughly 20 MeV [14], [15], for recent reports see [16], [17] . This value of T decreased substantially as compared to the quoted value one decade ago, T = 192 MeV [18], that used calculations with improved staggered fermions for various light to strange quark mass ratios in the range [0.05, 0.5], and with a strange quark mass fixed close to its physical value (although with an estimate for the string tension 10% larger than the usually quoted), while in [19] the crossover temperature was reported to be close to the present day value for the renormalized chiral susceptibility and about 25 [20,21].
The necessity to devise powerful measures of signatures for this crossover, that could also be useful in the experimental searches, has become a main objective. Besides the chiral condensate and its derivative with respect to the quark mass, the chiral susceptibility, used to probe the restoration of chiral symmetry, fluctuations and correlations in conserved charges have become tools to identify the transition from hadronic to quark-gluon degrees of freedom in the crossover region. In relativistic heavyion colliders experiments the ratios of such fluctuations are obtained in precision experimental studies for several collision energies as part of the RHIC beam energy scan program [22,23], which focuses on the search of the CEP. Chemical freeze-out parameters are then extracted within a canonical ensemble description of the data. These parameters lie below the freeze-out parameters extracted formerly from particle yields in the Hadron Resonance Gas (HRG) model [24] , [25,26] and the single freeze-out model [27][28][29]. A recent analysis using the HRG model to fit the net Kaon fluctuations at RHIC [30] provides experimental evidence that the freeze-out temperatures for strange hadrons could be 10 − 15 MeV higher than for the light ones.
Other measures resort to the equation of state (EoS) for determination of the trace of the energy momentum tensor and related specific heat and speed of sound. It has been pointed out a long time ago that the EoS near the QCD phase transition might be very soft as compared to an ideal pion gas [31][32][33] such that a "longest lived fireball" could be produced in relativistic heavy ion collisions [34] with ideal conditions to study signatures of the quark gluon plasma (QGP) as it goes through the stage of deconfinement. This softest point in the EoS is seen as a minimum in the velocity of sound, which measures the rate of change of the pressure with respect to the energy density. lQCD results show that at the minimum the energy density is only slightly above that of normal nuclear matter density [15].
In the present study we address the effects of the light and strange current quark masses in the light of an effective theory of QCD. As is well known from Chiral Perturbation Theory the canonical mass term represents only the leading order of an expansion in the masses themselves [35][36][37]. The explicit symmetry breaking pattern involves current quark mass dependent interactions at higher orders.
Our extension included first two kinds of eight quark chiral symmetry preserving interactions [53,54], which were needed to complete the number of vertices important for dynamical chiral symmetry breaking in four dimensions [55,56], and resolved instability issues related to the model's effective potential reported in [57]. The next extension added the set of explicit symmetry breaking (ESB) multiquark interactions at the specific N c order considered [58,59]. The phenomenological impact of the ESB terms on the quality of the low lying spectra of pseudoscalars and scalars as well as other related observables is remarkable, in comparison with the models without their inclusion. In particular the possibility to accurately describe the spectra of the scalar mesons together with a good fit for their strong decays gave us confidence that the model parameters obtained represented an adequate set for the analysis of the model related QCD phase diagram.
In a subsequent work [60] we have shown that the extended model leads to the emergence of two CEPs, associated with the light and strange quark condensates, in contrast to the common picture in which only the light condensate relates to a first order transition (except for a small effect on the strange condensate due to the coupling of both sectors), while the strange condensate displays a crossover behavior. The two CEPs act upon the onset for formation of strange quark matter, which is shifted to significantly smaller values of µ B .
Recently a similar extension of model interactions within the quark meson model and taking into account finite size effects led to the interesting result that quark matter with only u and d quarks may be the stable configuration in a region close to the end of the table of elements [61].
In face of these new developments one recognizes that the current quark mass effects are far from being fully explored and understood.
The crossover regime which follows after the two CEP's reported in [60] towards lower µ B until reaching the µ B = 0 line at higher temperatures is expected to weaken the prominent ESB features mentioned for the critical zone. Nevertheless we show that some observables are sensitive to the current quark mass values and ESB interactions. We stress that the numerical values of the current quark masses are intertwined with the dynamics of the ESB interaction terms and must always be considered together in the extended version. The numerical value of the strange current quark mass is reduced by roughly a factor 2 (thus being possible to reach its empirical value and light to strange quark mass ratios) when ESB interactions are considered in conjunction with best fits for the hadronic spectra, compared to the case without ESB interactions.
The obvious drawback of the model is the lack of confinement. However, it should be noted that the NJL model shares the global symmetries with QCD. Therefore it gives a reasonable tool to study the critical phenomena even if the location of the CEP may differ. For instance, both NJL and Polyakov-loop NJL models predict that the critical point is located inside the pion condensed phase. This result is consistent with the QCD no-go theorem [62], which is a rigorous statement in the large-N c limit. Thus, we suppose that symmetries combined with the 1/N c approach lead to a reasonable picture of the phase properties of the hadron matter even if the confinement mechanism is not included.
The comparison with lQCD data shows that there is room to improve the model calculations, mainly with respect to the transition to the Stefan-Boltzmann limit that occurs in a much shorter temperature interval, as well as a systematic shift to lower temperatures of the observables considered. The main observation is that the model is able to capture important features, such as the presence of a dip in the sound velocity, and the slopes of most observables considered.
Here we show that by coupling the model to the gluonic degrees of freedom through the Polyakov loop, the temperature gap between our model predictions and lQCD data is practically removed, as well as an improvement is obtained regarding the height of the peak of one of the correlators. An overall good agreement with lQCD data is obtained.
The text is organized as follows: after revisiting the model Lagrangian and thermodynamic potential in section II we present the thermodynamic observables that we compute and the model fits in section III,III A and III B and discuss the results in section IV. Our conclusions are summarized in section V.

A. Model thermodynamic potential
Although the model Lagrangian has been introduced and applied in previous works we indicate it here for completeness and refer for further details to [58][59][60]75].
The effective multiquark Lagrangian is expressed in terms of the U (3) Lie-algebra valued field Σ = (s a − ip a ) 1 2 λ a , involving the quark bilinears s a =qλ a q, p a = qλ a iγ 5 q; a = 0, 1, . . . , 8, λ 0 = 2/3 × 1, λ a are the standard SU (3) Gell-Mann matrices for 1 ≤ a ≤ 8. The q designates the color quark fields which enjoy the chiral flavor U (3) L × U (3) R global symmetry of QCD in the massless case. In addition the Lagrangian depends on external sources χ, which generate explicit symmetry breaking effects. In terms of these fields and sources the Lagrangian density reads to next to leading order (NLO) in the large N c counting with and the source dependent pieces Under chiral transformations one has for the quark fields R ; the sources transform as the field Σ. At this stage the sources can be fixed as the current quark masses χ = 1/2diag(m u , m d , m s ), after using the freedom related with the Kaplan-Manohar ambiguity associated with L 1 , L 9 , L 10 [76] and which are henceforth set to zero.
The couplings g i =ḡ i /Λ n (g i stands generically for any coupling) carry negative dimensions, given by the powers of Λ, and thus the Lagrangian is non renormalizable. Here Λ ∼ 4πf π ∼ 1 GeV [77] is associated with the scale of spontaneous chiral symmetry breaking.
As mentioned in the introduction this Lagrangian contains all non-derivative spin 0 multiquark interactions up to the same counting in large N c as the 't Hooft interaction, given by the term ∼ κ in Eq. 2. It has been shown that the large N c counting scheme selects the same interactions that are also relevant in the effective potential in the limit Λ → ∞, i.e. those scaling at most as Λ 0 . The LO contributions are only the 4q interaction ∼ G of the original NJL Lagrangian and the canonical quark mass term L 0 , all the other terms are N −1 c suppressed with respect to LO. These contain terms which violate the OZI rule (κ, κ 2 , g 1 , g 4 , g 7 , g 8 ), of which (κ, κ 2 ) break the U (1) A symmetry and are thus anomalous, as well as interactions which describe four-quark componentqqqq admixtures to theqq ones (g 2 , g 3 , g 5 , g 6 ).
The bosonization of the Lagrangian is carried over with functional integral techniques in the stationary phase approximation resulting in the effective mesonic Lagrangian density L bos at T = µ = 0, written in terms of the scalar, σ = λ a σ a , and pseudoscalar, φ = λ a φ a , nonet valued fields.
The result of the stationary phase integration at leading order, L st , is shown here as a series in growing powers of σ and φ. The coefficients h a , h ab , ... depend on the current quark masses and encode all the dependence in the coupling constants, see Eq. 12 below for h a . As in the case of the mass parameters, also only the h a with (a = 0, 3, 8) (or h i , (i = u, d, s) in the flavor basis) do not vanish [54].
The result of the remaining Gaussian integration over the quark fields is given by W ql . Here the second order operator in euclidean space is the constituent quark mass matrix resulting from the process of spontaneous symmetry breaking, which requires a redefinition of the field σ → σ + M, such that the new vacuum expectation value vanishes < σ >= 0. The one-quark-loop action W ql has been obtained by using a modified inverse mass expansion of the heat kernel associated with the given second order operator [78,79]. The procedure takes into account the differences ∆ ij = M 2 i − M 2 j , in a chiral invariant way at each order of the expansion, with b i being the generalized Seeley-DeWitt coefficients. The is the average over the regularized 1-loop euclidean mo-mentum integrals J i with i + 1 vertices (i = 0, 1, . . .) We use the Pauli-Villars regularization [80] with two subtractions in the integrand [81] We take only the dominant contributions to the heat kernel series, up to b 1 , b 2 for meson spectra and decays, which involve the logarithmically I 1 and quadratically I 0 divergent integrals in Λ.
The model thermodynamical potential Ω in the mean field approximation is written as a contribution stemming from the stationary phase approximation containing all the dependence on the model couplings, V st , and one which is related to the heat kernel quark one loop integrals J −1 which now carry the explicit T, µ dependence (for details please see [83]) where h i , i = (u, d, s) are solutions of the following system of cubic equations The h i are equal to one half the (unsubtracted) quark condensates, i.e. without the 2nd term in

III. THERMODYNAMICAL PROPERTIES OF STRONGLY INTERACTING MATTER
Strongly interacting matter is expected to undergo two transitions when subjected to high enough temperature (T ) and/or chemical potential (µ): deconfinement and chiral symmetry (partial) restoration. Although a straight connection between the two is still unclear they are for the most part expected to occur more or less simultaneously [84], [85].
The temperature and chemical potential dependence of fluctuations and correlations of conserved charges can be useful tools serving as indicators for the transition behavior. The fluctuations and correlations of the charges are given respectively by: Here we will consider the results pertaining to baryonic number N B , electric charge number N Q , and strangeness number N S . The corresponding chemical potentials are related to 1 As can be readily deduced from the relations of the corresponding The traced energy-momentum tensor Θ µ µ and the speed of sound C s are also thermodynamical quantities of interest. Both of these have been evaluated in lQCD thus we can use them as benchmarks to evaluate the adequacy of our models. They can be obtained respectively as: with P denoting the pressure, = T s−P the energy density, s = ∂P/∂T the entropy density and C V = (∂ /∂T ) V the specific heat at constant volume.

A. Polyakov loop extension
Although no gluonic degrees of freedom are present in the NJL model its extension to the so called Polyakov-Nambu-Jona-Lasinio Model is an attempt to mimic part of its dynamics by considering a static homogeneous background gluonic field in the temporal gauge The Polyakov loop L, winding around the imaginary time with periodic boundary conditions, and its trace in color space, φ (and charge conjugateφ) are given as where P stands for path-ordering and β = 1/T . In the quenched limit φ is an order parameter for the transition between the confined phase where the center of SU (N c ) symmetry (Z Nc ) is verified (vanishing traced Polyakov loop) and the deconfined phase where this symmetry is spontaneously broken [86]. An additional term, the Polyakov potential, must be added to drive this temperature induced spontaneous breaking. Its form can be determined by fitting lattice QCD observables. We take two choices, U I [66] and U II [67,74], with parameters shown in Table II, and N S = −Ns (by convention strangeness number is the negative of the number of strange quarks).
• Logarithmic form • Exponential K-Log form where the term proportional to K is the Van der Monde determinant, and the exponential term going with a 2 is a modification introduced in [74]. We considered a slight modification of this potential as we used U II φ,φ, T = U II φ,φ, T −U II [0, 0, T ] which enables the reproduction of the expected vanishing value for Ω/T 4 as we approach the vacuum ({T, µ} = {0, 0}). From a practical point of view the extension from NJL to PNJL amounts to the introduction of two new classical fields in the model, φ andφ, the introduction of the Polyakov potential and a modification of the occupation numbers (see for instance [72] for details on the implementation of the model) 2 . At vanishing baryonic chemical potential, the case considered in this work, φ =φ.

B. Parameter fitting
The parameters of the model are the quark current masses (m u , m d and m s ), the cutoff (Λ) and the couplings (G, κ, κ 2 , g 1 , g 2 , g 3 , g 4 , g 5 , g 6 , g 7 and g 8 ). In the formulation of the model without non-canonical explicit chiral symmetry breaking, NJLH8q, these 15 parameters are reduced to 8 (κ 2 and g i with i = 3, . . . , 8 are set to zero). If we neglect the isospin symmetry breaking (m l ≡ m u = m d ) these are further reduced to 7.
As was shown in [54] the NJLH8q model can be fitted for several fixed values of the OZI violating eightquark interaction coupling g 1 (the four-quark interaction strength, G, is smaller for increasing g 1 but the remaining parameters are unchanged) while keeping the mesonic spectra unchanged apart from a decrease in the σ meson mass for increasing value of g 1 (the scalar mixing angle also changes). Two sets, denoted as NJLH8qA and NJLH8qB, are shown in Table I with the latter corresponding to the highest value of g 1 (and conversely the lowest G). In these isospin symmetric sets (m l ≡ m u = m d ) we fit the model parameters by imposing a value of g 1 and fitting the remaining 6 parameters using the pion and kaon weak decay couplings (f π and f K ) and the meson masses (M π , M K , M η and M a0 ).
This freedom allowed us to isolate and study the impact of the eight-quark interaction term in the model phase diagram in the chemical potential-temperature, {µ, T }, plane [83]. One of the main highlights of this study was the realization that the CEP is shifted to lower chemical potential and higher temperature with increasing g 1 . This in turn leads to a substantial reduction in the related crossover temperature at µ = 0 compared to the case with weak g 1 coupling, as reported before in [6,89], with the lower values of T complying with lQCD results [19].
The extension to include the non-canonical explicit chiral-symmetry breaking interactions introduces 7 new parameters. For the parameter set NJLH8qmA from Table I we chose to impose the value of the current masses (m l and m s ). The remaining 12 parameters can be fitted by fixing f π , f K , the pseudoscalar and scalar mixing angles (θ ps and θ s ) and the 8 meson masses (M π , M K , M η , M η , M a0 , M K * , M σ and M f0 ), in the isospin limit. Note that the inclusion of the ESB interactions allows to fit the pseudoscalar as well as the scalar spectra to empirical data with a high degree of accuracy, as well as the weak decay constants and the current quark mass values. This on the other hand reduces the former freedom in the interplay of G, g 1 parameters, which is now considerably narrowed down, favoring the strong g 1 coupling strength, (we are considering here a range 470 MeV M σ 500 MeV). One also sees that the increase in g 1 comes in this case accompanied not only by a decrease in G, but also a decrease in the ESB couplings g 4 , g 7 .

IV. RESULTS AND DISCUSSION
A. Speed of sound and energy-momentum trace anomaly In Fig. 1(a) we see a comparison of the temperature dependence of the squared speed of sound at vanishing chemical potential as obtained using lQCD and several parametrizations of our model. Although a complete quantitative agreement seems impossible with the con- Table I. Model parameters obtained using a regularization kernel with two Pauli-Villars subtractions in the integrand (see [90] Parameters marked with an asterisk ( * ) were kept fixed. Several quantities which are either outputs or kept fixed (and used in the fit of the remaining parameters) are presented in the bottom rows: weak decay couplings ([fπ] = [fK ] = MeV), meson masses for the low-lying scalars/pseudo-scalars, the dynamical masses of the quarks (given in MeV) and the corresponding chiral condensates qq i (given in MeV 3 ) taking into account the subtraction of the contribution coming from the current mass, see Eq. 13. The pseudoscalar and scalar mixing angles (θps/θs) are given in degrees. Sets NJLH8qA and NJLH8qB include up to eight-quark interactions but no non-canonical explicit chiral symmetry breaking terms whereas set NJLH8qmA does include these terms.
Set m l ms G κ κ 2 g 1 g 2 g 3 g 4 g 5 g 6 g 7 g 8 Λ NJLH8qA 5.94 186. 12   sidered parametrizations an approximation of the general qualitative behavior can be obtained: the squared speed of sound increases with temperature starting at zero until a critical point upon which it dips into a local minimum (for 3 of the sets) and then rises again going asymptotically to the Stefan-Boltzmann limit. For temperatures above this local minimum the obtained speed of sound overshoots the lQCD result whereas for lower temperatures we obtain lower values. Note that as one expects the speed of sound to go to zero in the limit of vanishing temperature and chemical potential the lQCD results should be followed by a dip towards zero for lower temperatures.
The model velocity of sound displays a soft point close to the one of lQCD, at around T ∼ 150 MeV. For this feature to be present in the model, the parameter g 1 of OZI-violating 8-quark interactions is required to have a certain strength g 1 ∼ 3000 − 4000 GeV −8 , see set NJLH8qB without ESB interactions, and sets NJLH8qmA, NJLH8qmB with ESB in Table I; at the weak coupling g 1 = 500 GeV −8 of set NJLH8qA the relative minimum is absent and the velocity of sound shows a monotonous decrease. It had already been noticed some time ago that the strength of parameter g 1 has impact on the number of degrees of freedom [83]; a rather strong (more than 50 %) suppression of the artificial quark excitations (due to lack of confinement of the model) at T /T c > 0 was observed, in comparison to a PNJL model calculation [70]. Thus it is understandable that this property manifests itself in the occurrence of the relative minimum in the velocity of sound, as the model emulates partially the missing degrees of freedom attributed to the onset of deconfinement. The inclusion of the ESB breaking parameters does not change this important property, in spite of the fact that the accurate fit of the low lying spectra and related properties strongly constrains the parameters of the model in the vacuum. On the contrary: as mentioned before, the ESB interac-  Table I: dashed lines correspond to the NJLH8q sets whereas solid lines correspond to the NJLH8qm sets (respectively the ones without and with non-canonical explicit chiral symmetry breaking). The markers, labeled as WBEoS and HotQCDEoS, correspond to continuum extrapolated lQCD results taken respectively from [14] and [15].
tions together with the requirement of having good fits of the spectra rule out the smaller values of g 1 strengths. The fact that the former freedom in the model parameter g 1 is narrowed down, and specifically to the values that describe the soft region in the speed of sound, can be seen as a major result regarding the phenomenological importance of including ESB interactions in the model.
Apart from the relative minimum in the speed of sound, notice that the sets with ESB display a slight dip after the steep rise and before flattening out. This is a remnant of the two CEPs encountered in the model as described previously, and still visible at µ B = 0. This behavior may be guessed to be very subtly present in the lQCD points as well, although it would require further investigation to clarify this point.
Moving to Fig. 1(b) , one sees that the height of the peak of the trace of the energy momentum tensor is improved, as well as the slope before the transition, comparing with lQCD data, by the inclusion of the ESB terms (or the selection of the strong g 1 coupling without ESB). To understand that this is natural to expect within our model, we recall that the trace of the energy momentum tensor and the number of degrees of freedom ν(T ) = (90/π 2 ) P (T )−P (0) As mentioned above, the number of degrees of freedom for weak and strong g 1 coupling was obtained in [72,83], as function of T /T c , where T c is the crossover temperature. Converting by this factor one obtains the slope behavior displayed in Fig. 1(b). The slope after the peak is steeper in the model than in lattice calculations, but this is also expected, since the model approaches the Stefan-Boltzmann limit faster. Furthermore since in the non-Polyakov loop extended case there are no gluonic degrees of freedom present, this limit corresponds to a lower value (lim T →∞ −Ω N JL /T 4 = 31.5 π 2 /90 whereas lim T →∞ −Ω P N JL /T 4 = 47.5 π 2 /90 , see for instance [72,83]).
The energy density, , and pressure, P , as well as their derivatives with respect to the temperature, T , the specific heat, C V and entropy density, s, are depicted in Figs 2. Despite the reasonable agreement, apart from a shift towards lower temperatures, when we look at the energymomentum trace anomaly, Θ µ µ , (see Fig. 1(b)), in the individual thermodynamical quantities involved, P and , as well as s (See Figs. 2(a), 2(b) and 2(c)), the effect of the missing degrees of freedom is clearly present in their asymptotic behavior. In the cases with stronger eight-quark interactions (NJLH8qB, NJLH8qmA/B) the specific heat (see Fig. 2(d)) deviates from the lQCD with a marked peak around the transition region reflecting the faster transitional behavior. The slope of the curve for temperatures lower than the transition is however better reproduced in the cases with stronger g 1 .

B. Fluctuations and correlations of conserved charges
In Fig. 3 the fluctuations of several conserved charges are shown. The main gross feature for the baryonic susceptibility χ B 2 in Fig. 3(a) is that the slope improves significantly for the sets with strong g 1 coupling, in comparison with lQCD data. One sees further a noticeable change of slope between the steep rise and the flattening of the curves for sets including ESB interactions, NJLH8qmA and NJLH8qmB. A hint of such behavior using the parameter sets NJLH8q and NJLH8qm. The markers, labeled as WBEoS and HotQCDEoS, correspond to continuum extrapolated lQCD results taken respectively from [14] and [15].
seems to be present in the HotQCDEoS data as well. In a similar fashion do the slopes get improved in the fluctuations of strangeness χ S 2 , shown in Fig. 3(c) for the sets NJLH8qmA and NJLH8qmB. As opposed to these the slope is much steeper for these sets in the calculation of the electric charge fluctuations χ Q 2 displayed in Fig. 3(b), when compared to the lattice points. By expressing these fluctuations in terms of the quark number susceptibilities one sees that the weight of χ u 2 in χ Q 2 is larger than in χ B 2 by a factor 4. We have looked at the individual contributions within the model and found that the transition for χ u 2 occurs faster than for χ s 2 , as expected, and that the slope increases when the ESB interactions are taken into account, while the crossed contributions vanish, reflecting the fact that the model has no gluonic degrees of freedom [93,95], (the Polyakov loop introduces such a correlation, see section IV C below). So it seems that the slope of χ u 2 dominates the scene in χ Q 2 , due to the weighting factor, while the distribution of weights in the χ B 2 leads to the correct slope in comparison to lattice results. The χ S 2 provides for a clean probe of the strange quark number susceptibility, as this is the only contribution. For this case one sees that the lQCD slope is well reproduced with the ESB model sets.
Regarding the correlations displayed in Fig. 4, the same effect seems to be at work for correlation of baryonic and electric charges χ BQ 11 , shown in Fig. 4(a); it displays a too fast increase as compared to lQCD for the sets with ESB breaking. As opposed to this the correlations of baryonic and strangeness charges χ BS 11 in Fig. 4(b) show a slope in conformity with lQCD. The correlation of strangeness and electric charges also gets improved for . The markers, labeled as WBEoS and HotQCDEoS, correspond to continuum extrapolated lQCD results taken respectively from [93] and [94].  the ESB sets, Fig. 4(c). This can also be understood by looking at the dependence of these correlations on the quark number susceptibilities Since the correlations χ us 11 are smaller in magnitude in lQCD than the χ i 2 (i = u, s) [93] (χ ud 11 was shown to be small in the flavor SU (2) case [96]), and vanish identically in the model, the only dependence is on χ s 2 for χ BS 11 , χ QS 11 . We have seen that the slope for χ s 2 reproduces well the corresponding lQCD slope, which explains the satisfactory behavior of the slopes of χ BS 11 , χ QS 11 as well. The situation is different for χ BQ 11 , which depends on χ u 2 , for which a too steep slope occurred compared to lQCD.

C. PNJL extension
The impact of coupling the quark degrees of freedom to the gluonic sector, using the PNJL model extension with two types of potentials is discussed in this section.
The gross feature is a systematic shift of all the curves describing the observables of the last subsection to higher temperatures, which is an important effect in bringing most of the observables related with fluctuations and correlated charges closer to the lQCD curves.
However the effect on the velocity of sound and the trace of the energy momentum tensor depends strongly on the type of Polyakov loop potential used. Let us discuss first these observables.
In Fig. 5(a) the velocity of sound is displayed, calculated with the Polyakov loop potential U I in Eq. 17, showing that independently of the NJL parameter sets of Table I considered, a too deep relative minimum for the velocity of sound occurs, about a factor 2.5 smaller in magnitude, in comparison with lQCD. This result supersedes all the nuances discussed previously in relation with ESB terms. In Fig. 5(b) one sees that the peak of the trace of the energy momentum tensor is roughly twice the value of the lQCD one. These dominating characteristics are also present in the polylogarithmic variant of the Polyakov loop potential in [73], which are therefore not shown.
Contrary to this, the potential U II in Eq. 18 discriminates between the different NJL sets. Minima occur for the sets with ESB interactions (and large g 1 coupling without ESB) see Fig. 6(a), two shallow minima, and a more pronounced one for the ESB set with stronger g 1 coupling; the set PNJLH8qA without ESB terms corresponding to weak g 1 coupling does not display a minimum in the velocity of sound, as also verified in [97] using U II , and as it was the case without the Polyakov loop extension, see Fig. 1(a).
Regarding the trace of the energy momentum tensor it turns out to be fairly well represented in comparison with the lQCD results using U II , see Fig. 6(b). The individual thermodynamical quantities contributing to Fig.  6(b), and P , as well as their derivatives with respect to the temperature T , C V and s are depicted for the case of the potential U II in Figs 7. Overall, the correspondence of the presented quantities with lQCD is quite satisfactory (note that the inclusion of the extra degrees of freedom enables the correct asymptotic behavior for P , and s). The slight change of slope in around T = 0.18, compared to lQCD, results in a visible peak in C V . For C V and Θ µ µ the quark interaction parameter set which presents the best fit is NJLH8qA. For the other quantities the difference between quark interaction parameter sets (mainly in the transition region) is too small for the purpose of comparison with lQCD results. We omit the corresponding figures for the choice of U I since the calculations resulted in large deviations from the respective lQCD data, as one would expect from Fig. 5(b), and turn out not to be very instructive.
Turning to the fluctuations and correlations of the different charge numbers N B , N Q , N S one observes that all observables which had a good slope in the NJL model, i.e. χ B 2 , χ S 2 , χ BS 11 , χ QS 11 for the sets with ESB interactions, get shifted to higher temperatures and agree fairly well with the lQCD data for both Polyakov loop potential implementations, see Figs Unfortunately the slope of the correlator χ Q 2 is not improved with either potentials, see Figs. 8(b) and 9(b).
Finally we show in Fig. 12 that the coupling of the quark and gluonic degrees of freedom leads to a nonvanishing correlation between the light and strange quark numbers, albeit smaller than in lQCD, with the U II potential yielding a larger fraction. We also remark that this quantity is not sensitive to the details of the parametrizations in the quark sector.

V. CONCLUSIONS
We have used the three flavor NJL Lagrangian that has been enlarged in recent years to accommodate systematically current quark mass effects at NLO in the large N c counting scheme to address several thermodynamic observables. These explicit symmetry breaking (ESB) interaction terms are of the same order as the 't Hooft U A (1) breaking anomalous contribution and the previously introduced symmetry preserving eight quark interactions. It has been shown that the ESB terms play a very important role in the description of accurate characteristics of pseudoscalar and scalar mesons. This opened for the first time the possibility to study the model phase diagram of QCD with a set of parameters which reproduces the empirical spectra, together with current quark masses that fit the actual PDG values, allowing to narrow down the uncertainties related to the model parameters.
While the model reaches systematically the Stefan-Boltzmann limit too fast as compared to the lattice results and is systematically shifted to lower temperatures as compared to lQCD, there are some relevant features which are reproduced. We highlight the main results: (i) The ESB terms together with the realistic spectra select a region of parameters with strong OZI violating 8q coupling g 1 . We recall that without the ESB terms there was an interval of values for this coupling, which in an interplay with the 4q G coupling, left the spectra unchanged except for the σ(500) mass that got reduced for increasing g 1 . The freedom in g 1 was accompanied by a sliding CEP position in the model QCD phase diagram.
(ii) In the strong g 1 coupling regime enforced by the ESB terms the velocity of sound displays a soft point as predicted by relativistic heavy ion models and lQCD. This relative minimum is absent in the NJL model which contemplates only the 4q and 't Hooft interactions, or weak g 1 couplings.  Figure 5. In 5(a) the model Polyakov loop extension with potential UI in Eq. 17 is shown for the squared speed of sound as a function of temperature ([T ] = GeV) at vanishing chemical potential as obtained using the parameter sets from Table I. In the legend a "P" has been attached at the beginning of each parameter set, meaning that the Polyakov loop extension has been applied, PNJLH8qA and PNJLH8qB correspond to the sets without the ESB terms, PNJLH8qmA, PNJLH8qmB include the ESB interactions.The markers, labeled as WBEoS and HotQCDEoS, correspond to continuum extrapolated lQCD results taken respectively from [14] and [15]. In 5(b) is shown the energy-momentum trace anomaly for the same Polyakov loop potential.  Figure 6. In 6(a) and 6(b) the same observables are shown as in 5(a) and 5(b), but obtained with the Polyakov loop potential UII in Eq. 18. The markers, labeled as WBEoS and HotQCDEoS, correspond to continuum extrapolated lQCD results taken respectively from [14] and [15].
(iii) The trace of the energy momentum tensor displays a peak with height close to the lattice results; the slope of this quantity gets improved as compared to the model without the ESB interactions. However, although the strong g 1 coupling regime describes overall better slopes, it leads to a visible peak in the transition regime for the related quantity C V , which is not favored by lQCD data.
(iv) The slopes of the susceptibilities χ B 2 , χ Q 2 , χ S 2 of the conserved baryonic, electric and strange charges are sensitive to the weighting factors of the quark number susceptibilities χ u 2 , χ d 2 , χ s 2 that enter in their definition. We find that the slopes for χ B 2 and χ S 2 , as well as for the correlation involving these two charges, χ BS 11 , get substantially improved, while it is too steep for χ Q 2 . The observable χ S 2 is a clean probe for the slope of the strange quark susceptibility χ s 2 , which agrees well with the corresponding lQCD slope.
(v) Finally, by coupling the quark to the gluonic degrees of freedom via the Polyakov loop we observe that the temperature gap between the NJL and the lQCD curves disappears practically and the overall characteristics of the lQCD data is rather well . The markers, labeled as WBEoS and HotQCDEoS, correspond to continuum extrapolated lQCD results taken respectively from [14] and [15].
reproduced. For the trace of the energy momentum tensor the Polyakov loop potential U II is better suited to describe the lQCD data than the potential U I , within our model calculations. 2 ) and in 8(c) strangeness (χ S 2 ). Same notation is used as in 5(a) for the lines. The markers, labeled as WBEoS and HotQCDEoS, correspond to continuum extrapolated lQCD results taken respectively from [93] and [94]. ) and in 10(c) electric charge and strangeness (χ Q S 1 1 ). Same notation for lines as before. The markers correspond to continuum extrapolated lQCD results taken from [94] .  Figure 12. Correlation of the up and strange charges (χ u s 1 1 ) compared to continuum extrapolated lQCD results taken from [93] (labeled as WBEoS). As this quantity is almost insensitive to the parametrization of the quark interactions we chose to display the effect of the choice of Polyakov potential (for the fermionic part we chose set NJLH8qB). Without the inclusion of Polyakov loop this quantity vanishes.