Ju l 2 01 9 A hadron-quark hybrid model reliable for the EoS in μ B ≤ 400 MeV

We present a simple version of hadron-quark hybrid (HQH) model. The model is composed of the simple independent quark model for QGP states and an improved version of volume-excluded HRG model for hadronic states. The improved version of volume-excluded HRG model yields the pressure as a simple analytic form. The switching function from hadron states to QGP states in the present model has no chemical potential dependence. The present HQH model with the simple switching function is successful in reproducing the Polyakov loop at zero chemical potential and the EoS in μB ≤ 400 MeV.


I. INTRODUCTION
Lattice QCD (LQCD) provides a lot of information on hot QCD. In particular, the recent 2+1-flavor LQCD simulation [1] has confirmed that the chiral and the deconfinement transition are "crossover" at finite temperature (T ) and zero baryon chemical potential (µ B = 0), where the continuum and thermodynamic limits were carefully taken. In general, such crossover nature means that the chiral and the deconfinement transition are not completely decoupled and the transition temperatures depends on the choice of observables. In fact, observable-dependent transition temperatures T  [2][3][4][5][6][7][8][9]. The definition and the determination are essential for the investigation of the presence or absence of critical endpoint (CEP) in QCD phase diagram. Particularly in Ref. [9], the LQCD data disfavors the existence of the CEP in µ B /T ≤ 2 and T /T (∆ l,s ) c (µ B = 0) > 0. 9. Another important subject is to understand LQCD data on T dependence of the equation of state (EoS), especially for the crossover region (100 MeV < T < 400 MeV). In the region, hadrons are supposed to melt into the strongly-correlated quark-gluon plasma; however, the physical interpretation of such hadron-quark transition has not been established yet. The EoS including the hadron-quark transition is necessary for the analyses of relativistic nuclear collisions. For these reasons, a lot of QCD data have been accumulated [1][2][3][4][5][6][7][8][9][10][11][12][13].
As a complementary approach to LQCD simulations, we can consider effective models. This approach is useful for the prediction of the transition lines, the existence and the location of the CEP, and the EoS. In fact, a lot of predictions are made for these quantities. The hadron resonance gas (HRG) model is a simple model for hadronic matter. The hadron is treated as non-interacting gas, and all hadrons listed in Particle Data Book [14] are taken into account in the model. The HRG model remarkably reproduces LQCD data on the various thermodynamic quantities in T < 1.3T (∆ l,s ) c (µ B = 0) [11], which indicates that one cannot neglect the excited hadrons even above the chiral transition temperature and hence hadrons possibly coexistent with quarks and gluons. The simultaneous treatment of quarks and hadrons has been discussed for a long time. The various models have been proposed so far, such as the quark-meson model [15] and the Nambu-Jona-Lasinio (NJL) model with mesonic loops, quark-diquak coupling and chiral soliton, However, excited hadrons are absent in such models.
In our previous papers [20,21], we proposed the hadronquark hybrid (HQH) model in order to describe the coexistence of quarks and hadrons. In the model, the total entropy s(T, µ B ) is divided into hadron and quark pieces: Namely, where the function s H (s Q ) is the entropy density for hadronic matter (quark-gluon plasma). The weight function f H means the occupancy of hadronic matter in the total entropy and is constrained in 0 ≤ f H ≤ 1. The s(T, µ B ) was determined from LQCD data on T dependence of s LQCD and the secondorder susceptibilities at µ B = 0. We apply HRG model for s H and independent-quark (IQ) model for s Q . The IQ model is a simplified version of Polyakov-loop extended Nambu-Jona-Lainio (PNJL) model [22][23][24][25], that is, this model treats the coupling between the quark field and the homogeneous classical gauge field, but not the couplings between quarks. The neglect of quark-quark interactions are justified from the fact that the light-quark chiral condensate is quite small in T > 220 MeV where s Q > s H . In our previous version of HQH model [20,21], we have confirmed that the HQH model well describe the LQCD data on the EoS for both T ≤ T c and T ≥ T c [16][17][18][19][20][21]. As another advantage of our approach, s LQCD automati-cally satisfies the thermodynamic inequality and the Nernst's theorem [26], In the HRG model, the interactions between baryons (antibaryon) are neglected, but it should be taken into account for µ B dependence of thermodynamic quantities. A simple way of treating volume-exclusion effects (repulsive force) [27] was suggested in Refs. [28,29]. This model is called "excludedvolume HRG (EV-HRG) model". Furthermore, a method of treating an attractive force in addition to the repulsive force was proposed in Ref. [30]. The volume-exclusion effects are included by fitting the volume b = 4 · 4πr 3 /3 [26] to either LQCD data or the core radius r of nucleon-nucleon force [28,29]. In the framework of Refs. [28][29][30], the interaction between baryon and anti-baryon and the radius of meson are neglected.
In this paper, we improve the HQH model of Ref. [21], taking the EV-HRG model for the hadron piece and using the simple IQ model for the quark-gluon piece. The EV-HRG model taken yields the pressure as a simple analytic function and guarantees that the pressure is µ B even. We refer to the present version of HQH model as "simple HQH (sHQH) model".
The switching function f H is determined from s LQCD at µ B = 0, i.e., f H has no µ B dependence. The sHQH model with the switching function well accounts for the Polyakov loop at zero chemical potential and the EoS in µ B ≤ 400 MeV, without introducing µ B dependence to f H , where the core radius r = 0.335 fm is taken.
The ∆ l,s and the Φ signal the chiral and the deconfinement transition, respectively. The ∆ l,s calculated with the HRG model becomes negative in higher T [5], whereas the corresponding LQCD result is positive. The present model has the same problem. As an interesting result of LQCD simulations in Ref. [5], the peak position of d∆ l,s /dT agrees with that of dε/dT at µ B = 0. In Ref. [7], furthermore, the transition line is estimated by the peak of dε/dT for finite µ B . Therefore, we use the peak and the half-value width of dε/dT as a transition region in µ B -T plane. We guess that the transition region determined from ε is close to the chiral-transition region calculated with LQCD simulations [8]. As a result, we show that both the regions are close to each other in µ B ≤ 400 MeV.
As a deconfinement-transition region, we take the peak and the half-value width of dΦ/dT and predict the transition region for µ B ≤ 400 MeV. We also determine a transition line from isentropic trajectories, and show that the transition line is between the deconfinement line and the transition line determined from ε. This paper is organized as follows. In Sec. II, we show the model building. Numerical results are shown in Sec III. Section IV is devoted to a summary.

II. MODEL BUILDING
We improve the HQH model of Ref. [21], modifying the EV-HRG model for the hadron part and taking the IQ mode for the quark-gluon one.
For the 2+1 flavor system, we can consider the chemical potentials of u, d, s quarks by µ u , µ d , µ s , respectively. These potentials are related to the baryon-number (B) chemical potential µ B , the isospin (I) chemical potential µ I and the hypercharge (Y ) chemical potential µ Y as As for µ I and µ Y , the right-hand side of Eq. (3) comes from the diagonal elements of the matrix representation of Cartan algebra in SU (3) group:

A. HRG model
For later convenience, we start with the HRG model. In the model, the pressure P H is divided into the baryon (B) part P B , the anti-baryon (aB) part P aB and the meson (M) part P M : with is the mass and the chemical potential of the i-th baryon (j-th meson), respectively. Here we have used the shorthand notation for the inegration over 3d-momentum p. In Eq. (5), all the hadrons listed in the Particle Data Table [14] are taken.

B. Modified version of EV-HRG
We first explain the EV-HRG model of Refs. [28][29][30]. The pressure P EV;H is obtained by with Here the effective baryon and anti-baryon chemical potentials, µ EV:B,i and µ EV:aB,i , are defined by whereb = bT 3 for a positive volume parameter b. It is not easy to obtain P EV;B and P EV;B , since µ EV;B,i (µ EV;aB,i ) includes P EV;B (P EV;aB ). Actually, P EV;B and P EV;aB are obtained by solving Eqs. (11) and (12) numerically.
In QCD, the pressure is charge-conjugation even (µ B even). Hence the P EV;H should be µ B even, because it is a model of explaining QCD in T < T c . However, µ EV:B,i includes a µ Bodd term µ B and a µ B -even termbP EV;B /T 4 , so that P EV;H is not µ B even.
The P EV;B and P EV;aB are now modified as The sum of P mod;B and P mod;aB are µ B even, since the sum is invariant under µ B → −µ B . For this reason, we take Eqs. (15)- (18). These equations show that P B ≥ P aB . The P mod:B and P mod:aB can be rewritten into for LQCD data on the EoS are available for T ≤ 400 MeV and µ B ≤ 400 MeV [5,7]. We then consider this region. We consider P B , because of P B ≥ P aB . The ℓ convergence of Eq. (19) becomes worse as |(µ B − m i )/T | becomes larger; note that K 2 (x) is proportional to exp(−x) for large x and µ B − m i is negative. Therefore, the convergence is worst for the smallest case (939−400)/400 where T = µ B = 400 MeV and m N = 939 MeV. Taking the ℓ = 1 term only is a 3 % error in Eqs. (19). In actual calculations, nucleon contribution in P B is only 3 %, so that taking the ℓ = 1 term only corresponds to 0.1% error. We can identify P B with its ℓ = 1 term and P aB with its ℓ = 1 one. This approximation is called "ℓ = 1 identification in this paper Using the ℓ = 1 identification, we can rewrite P mod:B as Multiplying both the sides of Eq. (22) byb exp(bP mod;B /T 4 ) and using the ℓ = 1 identification, one can obtain Noting that the Lambert W (z) function is the inverse function of W e W = z, one can get P mod:B as a simple analytic function: Namely, In the limit ofb = 0, the P mod:B tends to P B , because of W (z) → z. Parallel discussion is possible for anti-baryon. The result is Hence the hadronic pressure becomes P mod;H = P mod;B + P mod;aB + P M with Eqs. (24) and (25). The entropy density s mod:H is obtained from P mod:H as This modified version of EV-HRG model is referred to as "modified EV-HRG model". Figure 1 shows T dependence of the total pressure P (T ) for µ B = 0, 400 MeV. The results of modified EV-HRG and HRG models are compared with LQCD ones [13]. In the modified EV-HRG model, we take the core radius 0.335 fm as a value of r, i.e., b = 0.64 fm 3 . For µ B = 400 MeV (lower panel), the EV-HRG result (solid line) agrees with LQCD one [13] in T < 210 MeV, while the HRG result (dashed line) is consistent with LQCD one in T < 150 MeV. For µ B = 0 MeV (upper panel), both the EV-HRG and the HRG result are consistent with LQCD one [13] in T < 210 MeV.

C. Independent quark model
We have to consider QGP states in the region T > 200 MeV by using the simple IQ model, since f H (T ) < 1, as shown later in Fig. 2. The Lagrangian density of the IQ model is where m f is the current mass of f quark and D µ = ∂ µ − igA a µ λa 2 δ µ0 with the Gell-Mann matrix λ a in color space. See Refs. [24,25] for the definition of the Polyakov loop Φ and its conjugateΦ.
Making the path integral over quark fields leads to where (29), the vacuum term has been omitted, since the pressure calculated with LQCD simulations does not include the term. The Φ andΦ are obtained by minimizing Ω Q = −P Q .
The entropy density s Q is obtained from P Q as We take the Polyakov-loop potential of Ref. [21]: The parameters a 0 , a 1 , a 2 , b 3 and T 0 were fitted to 2+1 flavor s LQCD in 400 < T < 500 MeV; see Fig. 1 of Ref. [21] for the fit. The resulting values are tabulated in Table I.
in the sHQH model, where it is assumed that the f H (T ) has no chemical-potential dependence. Note that s mod:H and s Q have chemical-potential dependence. Therefore, f H (T ) is determined so as to s = s LQCD [13] at µ B = 0: Namely, In Fig. 2, the f H (T ) of Eq. (37) is shown by dots with error bars. The errors come from s LQCD . The solid line is a fitting function for the f H (T ) of Eq. (37); in the χ 2 fitting, the line is assumed to be 1 in T < 180 MeV. From now on, we regard the solid line as the switching function f H (T ).
The pressure P with no vacuum contribution is obtainable from s LQCD of Eq. (36):

III. NUMERICAL RESULTS
As mentioned in Sec. I, we consider the transition region determined from with the peak and the half-value width of dε(T, µ B )/dT and the deconfiment-transition region with the peak and the half-value width of dΦ(T, µ B )/dT .
A. T dependence of the Polyakov loop for µB = 0 ∼ 400 MeV Figure 3 shows the Polyakov loop Φ as a function of T for the cases of µ B = 0, 100, 200, 300, 400 MeV. The LQCD result is available only for µ B = 0 MeV [5]. In the upper panel for µ B = 0 MeV, the sHQH result (solid line) well reproduces LQCD one in which the continuum limit is taken. We then predict the Φ for µ B = 100, 200, 300, 400 MeV in the lower panel. µ B dependence of Φ is small.

B. Transitions
We first consider the case of µ B = 0. Table II shows Figure 4 shows the transition region T ε c determined from the peak and the half-valued width of dε(T, µ B )/dT and the lattice chiral-transition region in µ B -T plane; the former is calculated with the sHQH model and the latter is analytic continuation of LQCD simulations from imaginary to real µ [8]. The transition region determined from dε(T, µ B )/dT is shown by a horizontal line with cross for each of µ B = 0, 100, 200, 300, 400 MeV; the cross is a maximum value of dε/dT and the line means the half-value width of dε/dT . The red solid line is made by connecting the crosses. Meanwhile, the blue band indicates the width of the chiral-transition region extrapolated from the imaginary-µ B region [8]. The model result is consistent with the LQCD result. Figure 5 shows In Fig. 6, the solid curve is a line connecting the points at which the curvature of isentropic trajectory becomes maximum. This figure suggests that a transition line can be estimated by n/s in µ B -T plane. Hence, the transition calculated with n/s may be deduced from relativistic nuclear collisions. There is no evidence of attractor of isentropic trajectory in the sHQH model.

C. The EoS
In Sec. III B, we considered the transition line determined from dε/dT , where ε(T, µ B ) = sT − P + µ B n.
One uses the Tolman-Oppenheimer-Volkof equations (TOV) equation to study structure of stars. The input EoS of the TOV equation is P and n. In addition, the isentropic trajectories, n/s=const, are important to study relativistic nuclear collisions. Therefore, we focus on P , s, ε, n.
In order to compare the present model with the previous model [21], we take the same assumption " f H (T ) has no µ B dependence", in the the previous model.   Figure 8 shows T dependence of s, P , ε, at µ B = 0 MeV. The solid and dashed lines are the results of sHQH and HRG-HQH models, respectively. Seeing s(T ), we find that the fitting of f H (T ) is good, since the sHQH result agrees with LQCD data [7]. Also for P and ε, the sHQH model agree with LQCD data. Comparing the results of sHQH and HRG-EV HQH models, we can find that EV effects are small for µ B = 0 MeV. Figure 9 shows T dependence of s, P , ε, n at µ B = 100 MeV. The solid and dashed lines stand for the results of sHQH and HRG-HQH models, respectively. The s, P , ε, n of sHQH model reproduce LQCD data [7]. Comparing the two lines, we ca see that EV effects are small still for µ B = 100 MeV.
Figures. 10-12 shows T dependence of s, P , ε, n for µ B = 200, 300, 400 MeV. The results of sHQH model well reproduces the LQCD data [7]. EV effects become large as µ B increases from 200 MeV to 400 MeV.  [7]; note that n is deduced from s, P , ε. LQCD sHQH HRG-HQH Fig. 11: T dependence of s, P , ε, n at µB = 300 MeV. See the the text for the definition of lines. LQCD are taken from Ref. [7]; note that n is deduced from s, P , ε.  [7]; note that n is deduced from s, P , ε.

IV. SUMMARY
We have improved the HQH model of Ref. [21], modifying the EV-HRG model [28,29] for the hadron piece and using the simple IQ model for the quark-gluon piece. The modified EV-HRG model yields the baryon and antibaryon pressures as simple analytic functions of Eqs. (24)- (25), and ensures that the pressure is µ B even.
We have determined the switching function f H from s LQCD at µ B = 0. The sHQH model with the switching function f H (T ) well accounts for LQCD data on the Polyakov loop at µ B = 0 MeV. This makes it possible to predict the Polyakov loop for µ B = 100, 200, 300, 400 MeV. The EoS calculated with the sHQH model is successful in reproducing the corresponding LQCD data in µ B ≤ 400 MeV, without introducing µ B dependence to f H , where the core radius r = 0.335 fm is taken. The switching function f H has also a simple form, since it has no µ B dependence.
As an interesting result of LQCD simulations for µ B = 0 [5], the peak position of d∆ l,s /dT agrees with that of de(T, µ B )/dT . In LQCD simulations for finite µ B [7], furthermore, a transition line is estimated by the peak of dε(T, µ B )/dT . We can then guess that the transition region determined from ε is close to the chiral-transition region calculated with LQCD simulations. In fact, we show that the transition region determined from dε(T, µ B )/dT almost agrees with the lattice result [8] on the chiral-transition region in µ B ≤ 400 MeV. This may make it possible to define a chiral-transition region in µ B -T plane with the peak and the half-value width of dε(T, µ B )/dT . As a deconfinement-transition region, we take the peak and the half-value width of dΦ(T, µ B )/dT . As for the deconfinement transition, we predict the transition region and confirm that the deconfinement-transition line is above the transition line determined from dε(T, µ B )/dT . In sHQH model, there is no evidence of attractor of isentropic trajectory. We have also found that the transition line determined from isentropic trajectories is between the deconfinement line and the transition line determined from dε(T, µ B )/dT . The transition determined from isentropic trajectories may be deduced from relativistic nuclear collisions.