Thermodynamics of strongly interacting matter in a hybrid model

The equation of state and fluctuations of conserved charges in a strongly interacting medium under equilibrium conditions form the baseline upon which various possible scenarios in relativistic heavy-ion collision experiments are built. Many of these quantities have been obtained in the lattice QCD framework with reliable continuum extrapolations. Recently the Polyakov$-$Nambu$-$Jona-Lasinio model has been reparametrized to some extent to reproduce quantitatively the lattice QCD equation of state at vanishing chemical potentials. The agreement was precise except at low temperatures, possibly due to inadequate representation of the hadronic degrees of freedom in the model. This disagreement was also observed for some of the fluctuation and correlations considered. Here we address this issue by introducing the effects of hadrons through the Hadron Resonance Gas model. The total thermodynamic potential is now a weighted sum of the thermodynamic potential of the Polyakov$-$Nambu$-$Jona-Lasinio model and that of the Hadron Resonance Gas model. We find that the equation of state and the fluctuations and correlations obtained in this hybrid model agrees satisfactorily with the lattice QCD data in the low temperature regime.


I. INTRODUCTION
The thermodynamic properties of strongly interacting matter under extreme conditions are being actively investigated both theoretically as well as experimentally. The lattice formulation of quantum chromodynamics (QCD) on discretized space-time provides a first principle approach in this direction [1]. In a system with light quarks a rapid crossover from color confined and chiral symmetry broken hadronic phase to a chirally restored and color deconfined partonic phase has been predicted [2][3][4][5][6][7][8][9][10][11][12]. In the physical case of two light and a heavier strange quark, lattice QCD simulations at zero net conserved charges, this cross-over temperature is expected to lie in the range 150 MeV < T c < 160 MeV as reported by the Hot-QCD [13,14] and Wuppertal-Budapest [15] collaborations. Though at a cross-over there is no singularity, it is observed that near T c various thermodynamic quantities exhibit a rapid change [16][17][18].
At the same time, it is also important to properly explore the regions of QCD phase diagram away from the temperature axis. It is expected that at some critical high baryon chemical potential and small temperatures the system may undergo a first order phase transition from hadronic to partonic phase [19][20][21]. This first-order phase boundary would continue for some lower chemical potentials and higher temperature, eventually terminating at a critical end-point (CEP) [22][23][24][25]. One of the major goals in the experiments with heavy ion collisions is to map this phase diagram of QCD and locate the CEP, if any.
A reliable way to understand the phase transition dynamics is through the study of correlations and fluctuations of conserved charges. At finite temperatures and chemical potentials fluctuations of conserved charges are sensitive indicators of the transition from hadronic matter to partonic matter. Moreover, the existence of the CEP can be signalled by the singularities in fluctuations. In lattice QCD framework many of these susceptibilities have been obtained at vanishing chemical potentials. Unfortunately for non-zero baryon chemical potentials the lattice QCD framework faces certain technical difficulties. Recently various intelligent techniques have been developed to circumvent these difficulties to some extent [4,7,8,[26][27][28][29][30][31][32][33][34][35][36][37].
A parallel approach with QCD inspired models are being developed alongside the lattice QCD approach to gain some insight into the various aspects of strongly interacting matter. Here we shall discuss one such model − the Polyakov loop enhanced Nambu−Jona-Lasinio (PNJL) model. Originally the PNJL model was introduced to enhance the Nambu−Jona-Lasinio (NJL) model [38][39][40][41][42][43][44] with the gluon thermodynamics effectively through the Polyakov loop. In the absence of the gluon dynamics, the NJL model does not incorporate the mechanism of confinement adequately. Extending the NJL model to the PNJL model by introducing a temporal background gluon-like field along with its self interactions, restores some sense of confinement into the model [45][46][47]. The NJL as well as the PNJL models conserve all the global charges like the chiral, baryon number, electric charge, strangeness etc., as in QCD. The multi-quark interactions in this model are responsible for the dynamical generation of mass, leading to spontaneous breaking of chiral symmetry.
In the initial parametrizations of the PNJL model the NJL parameters were set from the masses of observed hadron masses and decay coefficients, while the Polyakov loop parameters were set from the pure gauge dynamics on the lattice. This already led to qualitatively similar results as in lattice QCD framework [47][48][49][50][51][52][53][54]. Over the years several studies were done to analyze the properties of this model as well as to improve the model step by step. For example in order to stabilize the ground state of the 2+1 flavor system improvements were introduced in the PNJL model [55][56][57] following similar improvements in the NJL model [58][59][60][61]. A first case study of the phase diagram in β-equilibrium using the PNJL model was reported in [62]. In a related work [63], the SU(3) color singlet ensemble of a quark-gluon gas has been shown to exhibit a Z(3) symmetry and within stationary point approximation it becomes equivalent to the Polyakov loop ensemble. Significant effects on baryon-isospin correlations due to the mass difference between the two light quarks were reported in [64]. Similarly the effects of finite system sizes on various thermodynamic quantities including fluctuations and correlations of conserved charges have been reported in ref. [65,66]. Also the first model study of the net charge fluctuations in terms of D-measure from the PNJL model [67] has been reported. In fact the validity of the fluctuation-dissipation theorem has also been discussed in the context of the PNJL model [68]. Viscous effects may play important role in the evolution of the hot and dense system. Their effects in terms of transport coefficients have been done in the NJL and PNJL model [69][70][71][72][73][74][75][76] and compared with hadron resonance gas studies [77][78][79][80]. The behavioral pattern of various observables extracted from the PNJL model may be found in [62,[81][82][83][84][85][86][87][88][89][90][91][92], The QCD phase structure has also been investigated for imaginary chemical potentials in the PNJL model framework [93][94][95]. Different interesting features of the Polyakov loop have led to the development of different formalisms of the PNJL model [96][97][98][99]. Effects of consideration of gluon Polyakov loop have been discussed in ref. [100][101][102]. Recently improvements of the Polyakov loop potential have been carried out by introducing the effects of back-reaction of the quarks [103,104].
The qualitative agreement between the various observables computed in the PNJL model and that in lattice QCD has been quite satisfactory. This agreement seemed very convincing once the temperature dependent observables were plotted against T /T c , where T c in the model was not equal to that obtained on the lattice. However lattice QCD data used for these studies were only reported at finite lattice spacings. Recently continuum extrapolations for a number of observables have been reported from lattice simulations [13,14,[105][106][107][108]. Correspondingly the various thermodynamic properties of strongly interacting matter were investigated within the framework of the PNJL model by reparametrizing the Polyakov loop potential [109] . The focus was to ascertain a quantitative agreement of a variety of observables with the lattice data. The overall correspondence was satisfactory but not perfect. One of the regions of disagreement was in the low temperature region where hadronic degrees of freedom are dominant. This may be expected as these degrees of freedom are not adequately addressed in the PNJL model.
On the other hand, the hadron resonance gas (HRG) model [110] has been very successful in describing the hadron yields in central heavy ion collisions from AGS up to RHIC energies [111][112][113][114][115][116][117][118][119][120][121][122][123]. Along with that the susceptibilities of conserved charges calculated in lattice QCD have been well reproduced by HRG model [124][125][126][127] for temperatures up to 150 MeV. Also the region of large chemical potentials below the critical region can be studied using this model. Thus this model is quite suitable for describing the hadronic phase of strong interactions. The HRG model is based on the Dashen, Ma and Bernstein theorem [128] which shows that a dilute system of strongly interacting matter can be described by a gas of free resonances. Though the long range attractive part of the hadron interactions is taken care of by these resonances the short range repulsion are also important for the description of strongly interacting matter. Near the critical/cross-over region HRG calculations tend towards Hagedorn divergence which may be due to the absence of repulsive interaction. This repulsive part is incorporated through the excluded volume effects in the HRG and is commonly known as EVHRG model [126,[129][130][131][132][133][134][135][136][137][138]. EVHRG equation of states have also been used for the hydrodynamic models of nucleus-nucleus collision [139][140][141]. Recently fluctuations of conserved charges, using HRG model, have been studied in Ref. [142] including the experimental acceptances of pseudo-rapidity and transverse momentum. In Ref. [137], higher moments of net-proton multiplicity have been studied using EVHRG model and compared with STAR data. A similar study was done by some of us for net-proton, net-charge and net-kaon using the HRG and EVHRG models [127]. The effects of finite volumes on the hadron gas were discussed in [143] and effects of magnetic fields were presented in [144]. In fact as the higher moments are expected to be more sensitive to the phase transition, any deviation of experimental observation from the model results may be taken as an indication of new phenomena. Investigation in this direction was done with variation of collision centrality as well as collision energy in [145].
Unfortunately, there is no single model which describes both the low temperature and high temperature domain accurately. However, a possible mechanism to address this issue has been suggested in [146]. The idea is to make an interface between one model appropriate for the hadronic phase and another model appropriate for the partonic phase. For the hadronic sector the authors have used the HRG model, and for the partonic sector they have used the PQCD results, and have used a switching function at the interface. In a follow up work [147] they have discussed the various versions of HRG models along with the PQCD data and obtained quite a number of thermodynamic variables including some susceptibilities in the baryon sector. They found the equation of state to be reproduced very well vis-a-vis the lattice QCD data, but the values of some of the susceptibilities deviate slightly from the lattice QCD data. However the baryonic susceptibilities obtained in STAR data are well reproduced with the EVHRG model for the hadronic sector [148].
In a similar study an interacting quark (IQ) model, which is a modified PNJL model was used for the partonic phase in [149,150], and HRG model for the hadronic phase. Here the multiquark interactions in the NJL sector have been neglected and the quark masses are taken to be the nominal Lagrangian masses. In these studies the equation of state and susceptibilities of various conserved charges have been computed and reproduces the lattice QCD data quite well.
In the present study we shall explore this approach utilizing the PNJL model and the HRG model. Before going on to the details let us clarify some of our assumptions and some of the differences with the earlier studies. Here we shall consider the simplest non-interacting form of the HRG model to describe the low temperature phase. We expect the effects of the short range repulsive forces incorporated in the excluded volume HRG model to be small for T < 150 Mev [127]. As we shall see, the non-interacting HRG model easily describes most of the physics in the low temperature phase.
Secondly, in both sets of the earlier studies the switching function is considered to have both temperature and chemical potential dependence. This should be true in general. However for simplicity we have chosen it to be temperature dependent as in the present work we shall only concentrate on observables at zero chemical potential.
Lastly, all the earlier studies have considered derivatives of the switching function for computing further observables using the thermodynamic relations. For example if switching function is used as an interpolating function of hadronic and partonic pressure, the entropy will be obtained from the temperature derivative, which also acts on the switching function. We however assume the switching function has no thermodynamic origin. It is used as a global interpolator between the hadronic and partonic phases irrespective of the observables. So thermodynamic operators should not act on the switching function. Nonetheless, we have discussed what happens if the temperature derivatives of switching function are included for some of the thermodynamic observables. As we discuss below that with proper parametrizations the numerical results do not depend much on whether the derivatives are included or not. Similar derivatives of the switching functions were considered in the IQ-HRG approach [150], where the derivative terms were used as parameters to fit the lattice QCD data for entropy and susceptibilities. Since we shall not consider the chemical potential dependence of the switching function, these derivative effects in the susceptibilities in our study are zero. Yet we shall see that the results seem to satisfactorily reproduce lattice QCD data in the hadronic phase.
We have organized the manuscript as follows. In the next two sections II and III, we shall briefly discuss the PNJL model and the HRG model respectively. In section IV we discuss how to couple the PNJL and HRG models, utilizing some of the thermodynamic quantities relating to the equation of state. Thereafter in section V we discuss the predicted behavior of the fluctuations and correlations of conserved charges from the PNJL-HRG model. And finally in section VI we summarize and conclude.

II. PNJL MODEL
In this section we briefly discuss the PNJL model used in this study. A detailed discussion may be found in Ref. [109]. The PNJL model was initialized with a Polyakov loop effective potential being added to the NJL model [45][46][47]. While the chiral properties are taken care of by the NJL part, the Polyakov loop explains the deconfinement physics. Extensive studies have been carried out using PNJL model with 2 and 2+1 flavors [47-49, 51, 55, 151-155].
Here we consider 2+1 flavor PNJL model taking up to six and eight quark interaction terms as in [55,151]. The thermodynamic potential in this case reads as, The fields σ f = ψ f ψ f correspond to the two light flavor (f = u, d) condensates and the strange (f = s) quark condensate respectively. The model has a four quark coupling term with coefficient g S , a six quark coupling term breaking the axial U(1) symmetry explicitly with a coefficient g D , and eight quark coupling terms with coefficients g 1 and g 2 necessary to sustain a stable minima in the NJL Lagrangian. The corresponding quasiparticle energy for a given flavor f is E f = p 2 + M 2 f , with the dynamically generated constituent quark masses given by, In the above, if σ f = σ u , then σ f +1 = σ d and σ f +2 = σ s , and so on in a clockwise manner. The finite range integral gives the zero point energy. The different parameters as obtained from [55] are given in Table I.  The finite temperature and chemical potential contributions of the constituent quarks are given by the next two terms. Note that these are basically coming from the fermion determinant in the NJL model modified due to the presence of the fields corresponding to the traces of Polyakov loop and its conjugate given by Φ = T rcL Nc andΦ = T rcL † Nc respectively. Here L( x) = Pexp i is the Polyakov loop, and A 4 is the temporal component of background gluon field.
The effective potential that describes the self interaction of the Φ andΦ fields are given by U ′ . Various forms of the potential exist in the literature (see e.g. [50,51,83,156,157]). We shall use the form prescribed in [51] which reads as, Here U(Φ,Φ, T ) is a Landau-Ginzburg type potential commensurate with the global Z(3) symmetry of the Polyakov loop [47]. J(Φ,Φ) is the Jacobian of transformation from the Polyakov loop to its traces, and κ is a dimensionless parameter which is determined phenomenologically. The effective potential is chosen to be of the form, The coefficient b 2 (T ) is chosen to have a temperature dependence of the form,  and b 3 and b 4 are chosen to be constants. The set of parameters is given in Table II.
Using realistic quark masses the deconfinement temperature obtained in lattice QCD is found to be much higher than the chiral transition temperature [9,11]. However the deconfinement temperature as measured from the peak of the entropy of a static quark is found to be consistent with the chiral transition temperature [158]. In our model framework we consider the temperature derivatives of the mean fields and locate their peaks to obtain the transition temperature. The temperature derivative of the Polyakov loop is closely related to the definition of temperature derivative of the static quark free energy that gives its entropy as defined in [158]. The corresponding T c was obtained from the average of the two peak positions. The resulting values of T c are listed in Table III.
With this parametrization we have been able to achieve both a crossover temperature of T c ∼ 160 MeV as well as quantitative agreement of temperature variation of pressure and various other observables with the lattice QCD continuum estimation [109]. However the quantitative agreement though close, was not precise enough. Various observables showed disagreement in different ranges of temperatures. The most common discrepancies were found in the low temperature region where the hadronic degrees of freedom dominate. We shall thus make an attempt to remove this lacunae by coupling the PNJL model with the HRG model.

III. HADRON RESONANCE GAS MODEL
We now discuss the HRG model briefly. The detailed discussions may be found in Refs. [110-115, 117, 118, 120, 126, 129-133]. The grand canonical partition function of a hadron resonance gas [110,126] can be written as, where the sum is over all the hadrons. id refers to an ideal gas of hadronic resonances. The partition function for the i th resonance is, where V is the volume of the system, g i is the degeneracy factor, T is the temperature, E i = p 2 + m 2 i is the single particle energy, m i is the mass and are respectively the baryon number, strangeness and charge of the particle, µ , s being corresponding chemical potentials. The (+) and (−) sign corresponds to fermions and bosons respectively. From the partition function one can calculate various thermodynamic observables of the system created in heavy ion collisions. In this work we shall incorporate all the hadrons listed by the Particle Data Group [159] up to mass of 3 GeV.

IV. COUPLING THE HADRONIC AND PARTONIC SECTORS
Let us now discuss the framework in which the hadronic matter described using the HRG model can be smoothly switched to partonic matter modelled through the PNJL one. The basic procedure is similar to the ones reported in [146][147][148][149][150]. The pressure of the system as taken to be a sum of partial pressures of the hadronic and partonic matter, weighted with a switching function as, where P P (T ) and P H (T ) are the pressures of partonic and hadronic sectors respectively and S(T ) is the switching function. For the partonic pressure we shall make a comparative study between PNJL models with six and eight quark type interactions respectively. Ideally S(T ) should be zero in the hadronic phase and 1 in the partonic phase. However since there is no singular phase boundary separating the two phases, we consider the switching function to interpolate smoothly from 0 at low temperatures to 1 at high temperatures. We have assumed the switching function to be independent of chemical potential at zero chemical potentials. The functional form of the switching function is given as, Here T S and ∆T S (T ) are parameters whose values should be closely related to the crossover temperature and width of the crossover region respectively. Note that the form of the switching function is somewhat different from the earlier studies. The essential difference is the introduction of two temperature scales − one being the central value T S of the switching and the other being the spread ∆T S . T S is related to the cross-over temperature and ∆T S (T ) to the spread of the switching function on the temperature axis from its central value. Given that the cross-over temperature is not sharply defined, we have let T S to be fixed at 160 MeV. Thereafter we try to find the best possible value of ∆T S (T ) so as to obtain the closest quantitative agreement of pressure as obtained in lattice QCD. The procedure for obtaining the various thermodynamic observables are as follows. We first obtain the pressure in the HRG model as well as the PNJL model from the respective thermodynamic potentials. Thereafter we obtain the hybrid pressure according to Eq. (8). Once we have the pressure the other thermodynamic quantities like entropy, specific heat etc., are obtained from the temperature derivatives of the hybrid pressure. In all the previous studies the temperature and chemical potential derivatives of the switching function in the hybrid pressure were included as well. We however consider the switching function to be defined globally for all the observables considered, so that the effective function should be independent of the observable being measured. But we shall still make a comparative study by both including and excluding the temperature derivative terms for some of the thermodynamic variables.
The entropy density is obtained from the first derivative of pressure with respect to temperature and is given by, Here s P (T ) and s H (T ) are the entropy density for the PNJL model and HRG model respectively. Note that the term in the square brackets denote the hybrid entropy due to the switching function, while the other term denotes the effects of the derivatives of the switching function itself. Similar would be the case for any other thermodynamic quantity where temperature derivative is involved. For example the energy density is given by, Here again the term within the square brackets would give the hybrid energy density and the other term comes from the derivative of S(T ). The derivative of S(T ) is given as, The derivatives of S(T ) contribute more terms for higher order temperature derivatives. For the specific heat at constant volume we get the expression as, Here the second order derivative term of S(T ) is given as,  We now show the temperature variations of some of the thermodynamic quantities. We shall compare the results in PNJL model with six and eight quark interactions respectively. We shall also compare results including and excluding the temperature derivatives of the switching function. As mentioned earlier the only parameter adjusted in the switching function is ∆T S (T ). It is allowed to take two different values on the two sides of T S when the derivatives of S(T ) are excluded. Here the values of ∆T S (T ) are not uniquely determined by any standard fitting procedure. However for the purpose of this paper it suffices to choose the values that best describe most of the thermodynamic quantities in the hybrid systems by trial and error. These values of ∆T S (T ) are given in Table (IV). Results are not much affected by shifting these values by 10%. The form of the switching function is shown in Fig. (1).  In Fig. (2) we show scaled pressure as a function of temperature for various models and compare them to the lattice QCD data [14,15]. While the HRG model agrees with lattice QCD data below T = 200 MeV, the PNJL model agrees with lattice QCD data above T = 150 MeV. On the other hand the combined PNJL-HRG model now agrees with the lattice QCD data for the full temperature range. The differences in the results due to the type of PNJL model chosen are quite insignificant. However some of those effects are absorbed in the values of ∆T S as given in table (IV). We note here that the pressure itself only depends on S(T ) and not on any of its derivatives. However the other thermodynamic variables are obtained from the derivatives of the hybrid pressure. As mentioned earlier, we are choosing ∆T S such that the bulk thermodynamic quantities obtained in lattice QCD framework are well reproduced in the PNJL-HRG model. Thus the values of pressure would depend on which form of S(T ) we are choosing. This is why we also show the plot for pressure for those values of ∆T S that are chosen when derivatives of S(T ) are included in computing rest of the thermodynamic variables.  The two other quantities that we considered for critically assessing the quantitative agreement with lattice QCD data are the trace of the energy-momentum tensor and the squared speed of sound. The trace of the energy-momentum tensor is defined as Θ µµ = (ǫ − 3P ). It is expected to be zero in a conformal theory. However the introduction of any scale like the cross-over or transition temperatures in the theory due to quantum interactions would break conformality and the trace will be non-zero. This is expected to be a very sensitive measure, inherent to the theory itself. Therefore a satisfactory agreement of the results obtained from lattice QCD with any other model is imperative for validation of the model as a simpler version of the theory.
In Fig. (3) we show scaled Θ µµ as a function of temperature for various models and compare them to the lattice QCD data [14,15]. Here the HRG model agrees with lattice QCD data up to about T = 170 MeV and rises much faster for higher temperatures. The PNJL model agrees with lattice QCD data above T = 150 MeV, except near the peak at around T = 200 MeV. Apart from this peak region, the combined PNJL-HRG model now agrees with the lattice QCD data for the full temperature range. The differences in the results due to the type of PNJL model chosen is again quite small. However near the peak region the eight-quark version of the PNJL model better describes the lattice QCD data and so does the PNJL-HRG model. The deviation from the lattice QCD data near the peak region would need further improvement of the PNJL model as the switching function favors it over the HRG model in this temperature range. Since our motivation here is to improve the situation for the temperature region dominated by hadrons, we do not intend to address this discrepancy here.
We now discuss the behavior of the speed of sound as a function of temperature. The squared speed of sound is given by, The variation of the squared speed of sound with temperature is shown in Fig. (4) for various models and compared with the lattice QCD data [14,15]. As seen from Eq. 14, the speed of sound on one hand contains the information of the equation of state and on the other hand contains the information of entropy to specific heat ratio. Thus the characteristics of the speed of sound depends strongly on the phase of strongly interacting matter. One expects that at very low temperatures the speed of sound would be small as the pressure of the hadronic system with masses much larger than the system is negligible. With increase in temperature the speed of sound will increase. However with increasing temperature the hadron resonances with higher and higher masses would be excited and the speed of sound may not reach the SB limit. In fact it may even start decreasing with temperature [126]. Using the HRG model all these features are obtained as shown in Fig. 4. After the crossover, the degrees of freedom are supposed to change from hadronic to partonic ones and therefore speed of sound may again increase. The minimum of the speed of sound known as the softest point may be a crucial indicator of the transition to be observed in heavy-ion collisions [160]. Such a minimum in the temperature variation of speed of sound is visible in the lattice QCD data as shown in Fig. 4, but is clearly absent in the PNJL model results. It is interesting to see that even the HRG model has a soft point of the same order of magnitude but at a higher temperature compared to the lattice QCD data. While the HRG model agrees with lattice QCD data up to about T = 150 MeV, the PNJL model results are consistent with the lattice data above T c . Thus the combined PNJL-HRG model shows reasonable agreement with the lattice QCD data in the full range of temperatures.
From the above discussions we may conclude that a hybrid model like the PNJL-HRG satisfactorily describes the equation of state of strongly interacting matter in a wide range of temperatures. We would now like to investigate whether the same holds true for the various fluctuations and correlations at vanishing chemical potentials. In strong interactions the net number of various quark flavors are conserved. But since quark states are not physically observed, one has to relate these flavor conservations into conservations of the experimentally observed charges of hadrons like baryon number B, electric charge Q and strangeness S. The fluctuations and correlations of conserved charges depend significantly on the state of strongly interacting matter at high temperatures and densities [17,56,57,161,162]. Variations of these quantities with temperature and various chemical potentials are supposed to carry signatures of phase transition or crossover [16, 65-67, 127, 143, 163-167]. They are related to the susceptibilities of pressure via fluctuation dissipation theorem [68]. The susceptibilities at various orders are easily obtained as a Taylor series expansion of pressure in terms of the corresponding chemical potentials. The diagonal Taylor coefficients c X n (T ) (X = B, Q, S) of n th order for the scaled pressure P (T, µ B , µ Q , µ S )/T 4 may be written in terms of the fluctuations χ X n (T ) of the corresponding order as, where the expansion is carried out around µ B = µ Q = µ S = 0. Similarly, the off-diagonal coefficients c X,Y n,m (T ) (X, Y = B, Q, S; X = Y ) of the (m + n) th order in the Taylor expansion of scaled pressure are related to the correlations between the conserved charges χ X,Y n,m (T ) as, At zero chemical potentials some of these fluctuations and correlations have been measured in the lattice QCD framework either in the continuum limit [12,105,106,[168][169][170][171] or for small lattice spacings close to the continuum limit [107]. We shall now discuss how much these quantities obtained in the PNJL-HRG model agree with the lattice QCD data. In the PNJL model the fluctuations and correlations are obtained by a suitable Taylor series fitting as has been discussed in detail in [48]. On the other hand these quantities are obtained easily in the HRG model by taking the corresponding derivatives with respect to the chemical potentials. Once the fluctuations and correlations are obtained in both the models, they are combined into the PNJL-HRG model using the same switching function S(T ) identified in last section. They are given by, and c X,Y m,n = S(T )c X,Y m,n P + (1 − S(T ))c X,Y m,n H As mentioned earlier we have assumed the switching function to be independent of the chemical potential in the close vicinity of zero chemical potentials. Let us now discuss some of these fluctuations and correlations obtained in the models and compare them to the lattice QCD data. In Fig. 5 the variation of the baryon number susceptibility c B 2 is shown as a function of temperature for various models and compared with lattice QCD data. Here we find c B 2 obtained in HRG model agrees with lattice QCD data up to about T ∼ 170 MeV. The PNJL model also agrees with lattice QCD quite well in this temperature range and therefore also the PNJL-HRG model. At higher temperatures the hybrid model is dominated by the PNJL model dynamics. There is a slight overestimation in the PNJL model which is also reflected in the PNJL-HRG model. Apart from that, the results in the six-quark and eight-quark versions of the PNJL models are numerically commensurate. The baryonic contribution of the PNJL model therefore seems to be sufficient in describing the strongly interacting matter even in the low temperature region described either by HRG model or the lattice QCD data. The variation of the electric charge susceptibility with temperature is shown in Fig. 6. There is a significant difference between PNJL model and lattice QCD results for c Q 2 below the crossover temperature T c . The lattice data is much larger than the PNJL model results. Though the baryon fluctuations in the lattice data are well accounted for by the constituent quarks in the PNJL model, proper considerations of other hadronic degrees of freedom below T c is crucial to obtain the correct values of electric charge fluctuations. On the other hand lattice QCD data is reproduced by the HRG model very well for T < 150 MeV. This is expected as the charge sector has dominant contributors from the light hadrons, which are practically absent in the PNJL model. Once the PNJL and HRG models are combined using the switching function, the PNJL-HRG model results again agree with the lattice QCD data very well.
The temperature variation of the strangeness susceptibility c S 2 is shown in Fig. 7. The computations in the HRG model seem to agree with the lattice QCD data up to much higher temperatures T ∼ 190 MeV. This is surprising given that the crossover temperature, as well as the temperature around which almost all other quantities computed in the HRG model start deviating from lattice QCD data at around T ∼ 160 MeV. One would expect the c S 2 in the HRG model to rise much faster and start deviating from lattice QCD data at much lower temperatures. One of the possible reasons for such a result may be that HRG model is constructed from experimentally observed hadrons, whereas the lattice QCD formulation could have contributions from additional species of strange hadrons, which are also predicted by quark model calculations [108]. At the same time the quantitative results for c S 2 obtained in PNJL model are significantly different from lattice QCD data up to T ∼ 250 MeV. The constituent masses of the strange quarks in the PNJL model have values above 500 MeV for T < 150 MeV. This is consistent with the omega baryon mass. However as the temperature rises, the constituent masses of strange quarks do not fall as fast as that of the light quarks. Therefore the strangeness fluctuations on their part, do not rise as fast as that of the light quarks as observed in the lattice QCD data. This aspect of the PNJL model would need further scrutiny and will be discussed elsewhere. Note that the range of temperature where the PNJL model disagrees with lattice QCD data is above the range where the switching from HRG model to PNJL model takes place. Inclusion of HRG model is therefore insufficient to address this issue. As a result the hybrid PNJL-HRG model does not agree well with the lattice QCD data. We thus encounter the first observable where the PNJL-HRG model could not satisfactorily describe the lattice QCD data in the full temperature range. It is apparent that this may be the case for other quantities that are strongly dependent on the strangeness content of the PNJL model. We now discuss the leading order correlations between the conserved charges. The correlator c BQ 11 between the baryon number and electric charge is shown in Fig. 8. In the hadronic phase the baryon and electric charges remain correlated as the baryons have positive electric charge and anti-baryons have negative electric charge. For small temperatures the correlations are small due to the relatively large masses. With increasing temperature however the correlation becomes non-zero. On the other hand for the 2+1 flavor system there are three quarks with equal baryon number in the partonic phase, but electric charge of down and strange quarks are together opposite of that of the up quark. At large temperatures when the quark masses are small with respect to the temperature, the BQ correlation should again tend to zero. In HRG model the BQ correlations keep on increasing as higher and higher mass states are getting excited with increase in temperature. On the other hand for lattice QCD data and PNJL model as well as any other model having a transition from hadronic to partonic phases, c BQ 11 would show a hump around the crossover region. This is shown in Fig. 8. Here the BQ correlation in the PNJL model is larger than that obtained in the lattice QCD data for T > 150 MeV. It may again be anticipated that this is due to the slow decrease of the strange quark mass with temperature in the PNJL model. The number of strange quarks is therefore much smaller than that of the light quarks, and they cannot compensate the electric charge of the up quarks sufficiently. The baryon number to strangeness (BS) correlation c BS 11 is shown in Fig. (9) and the electric charge to strangeness (QS) correlation c QS 11 is shown in Fig. (10). For the intermediate temperatures the BS correlation is a little less than the QS correlation due to the contributions from lighter strange mesons in the latter. In the HRG model, these two correlations again keep on increasing indefinitely as higher and higher hadronic states are getting excited. The behavior in the lattice QCD data as well as the PNJL model is as it should be in a theory that have partons in the high temperature phase. This can be understood by noting that at low temperatures the correlators are small due to the large hadronic masses. They will increase with temperature, and in the partonic phase the correlation saturates as there is only one species of quarks containing both strangeness and baryon number or strangeness and electric charge.
As in the case of c S 2 , the BS and QS correlators in the PNJL-HRG model underestimate the lattice QCD data in the intermediate range of temperatures possibly due to the slowly decreasing strange quark mass in the PNJL model. Therefore it seems that this deviation would pervade all other observables related to strangeness. A reparametrization of the NJL part of the model may be able to address this issue and will be discussed elsewhere.

VI. CONCLUSION
In this study we discussed a scheme to address the inadequacy of the PNJL model in describing the hadronic state of matter as pointed out in our earlier work in [109]. A straightforward approach would be to add the hadronic contribution from the HRG model. However a simple addition of HRG model to PNJL model would have led to overcounting the degrees of freedom. An interpolating function was therefore used in line with some earlier studies [146][147][148][149][150], to smoothly switch between the hadronic and partonic matter. In general this function would be dependent on temperature and chemical potentials. Since we have discussed various thermodynamic observables vis-a-vis lattice QCD data at zero chemical potentials, we have considered only a temperature dependent switching function.
In the earlier studies with the switching function, its temperature and chemical potential derivatives were considered for computing observables from thermodynamic relations. However we considered the switching function to be independent of the observable chosen. Otherwise the various derivatives of the switching function may be obtained by fitting with as many number of observables from lattice QCD see e.g. [150]. This would increase the number of free parameters in the model, thereby reducing its predictive power. Instead we emphasized on the better agreement of the PNJL model and lattice QCD data above the crossover region. Therefore a considerable part of the equation of state in continuum lattice QCD data, both below and above the crossover temperature, was well described by HRG and PNJL models respectively. Hence we only needed to interpolate in the crossover region using the switching function. Thereafter all other observables were obtained from thermodynamic relations acting on the individual pressure of the two models weighed with the switching function. For completeness however we have made a comparative study of observables related to the equation of state by both including and excluding these derivatives. All the thermodynamic quantities computed in the PNJL-HRG model reproduced the lattice QCD data quite well.
Once the switching function was fixed, we computed the predicted behavior of various fluctuations and correlations of conserved charges. We found that the fluctuations of baryon number and electric charge computed in the PNJL-HRG model is in good quantitative agreement with the lattice QCD data. However the strangeness fluctuations in the PNJL model was somewhat different from the lattice QCD data around the crossover region. As a result the PNJL-HRG model did not reproduce lattice QCD data satisfactorily. We emphasize that this disagreement is not due to the switching function but rather the inadequacy in the PNJL model. We argued that since in the PNJL model the constituent mass of the strange quarks does not decrease fast enough with the rising temperature, there is a departure of fluctuations of strangeness from the lattice QCD data. It naturally followed that none of the leading order correlators in the PNJL-HRG model could reproduce lattice QCD data as well.
Finally we would like to conclude that the scheme of introducing the HRG model to the PNJL model using the switching function seems to have worked well. Discrepancies if any are only due to difference between PNJL model and lattice QCD data in the relevant temperature ranges. Further for treading into the non-zero chemical potential regions we may have to consider a chemical potential dependence of the switching function. We hope to address these issues elsewhere.