Thermodynamic potential for quark-qluon plasma with finite quark mass and chemical potential

We summarize the derivation of the finite temperature, finite chemical potential thermodynamic potential in the bag-model approximation to quantum chromodynamics (QCD) that includes a finite s-quark mass in the Feynman diagram contributions for both zero-order and two-loop corrections to the quark interaction. The thermodynamic potential for quarks in QCD is a desired ingredient for computations of the equation of state in the early universe, supernovae, neutron stars, and heavy-ion collisions. The 2-loop contributions are normally divergent and become even more difficult in the limit of finite quark masses and finite chemical potential. We describe various means to interpolate between the low and high chemical potential limits. Although physically well motivated, we show that the infinite series Padé rational polynomial interpolation scheme introduces spurious poles. Nevertheless, we show that lower order interpolation schemes such as polynomial interpolation reproduce the Padé result without the presence of spurious poles. We propose that in this way one can determine the equation of state for the two-loop corrections for arbitrary chemical potential, temperature and quark mass. This provides a new realistic bag-model treatment of the QCD equation of state. We compute the QCD phase diagram with up to the two-loop corrections. We show that the two-loop corrections decrease the pressure of the quark-gluon plasma and therefore increase the critical temperature and chemical potential of the phase transition. We also show, however, that the correction for finite s-quark mass in the two-loop correction serves to decrease the critical temperature for the quark-hadron phase transition in the early universe.


Introduction
A description of the equation of state (EoS) of matter formed in the early universe, heavy-ion nuclear collisions, supernovae, or neutron stars should include the consequences of a possible phase transition between hadronic matter and quark gluon plasma (QGP). For many astrophysical applications the description of quark matter has only been considered in the zero-order bag model (e.g. [1][2][3][4][5][6]). However, the 2-loop corrections at the next order are also significant and should be considered in the EoS [7]. In this approach one can then construct the QGP EoS from a phase-space integral representation over the scattering amplitudes. Another approach is to model the QGP in the context of a chiral effective field theory (e.g. [8,9] and refs. therein). Although considerable recent progress has been made in this later approach, here we consider the former bag model approach because it is often employed in astrophysical simulations for its simplicity.
In this approach, it is convenient to compute the EoS for the QGP in terms of the grand thermodynamic potential, Ω(T, V, μ) [10]. Adopting the convention of Landau and Lifshitz [11], the thermodynamic potential can be defined in terms of the partition function Z as: where b º kT 1 (henceforth we adopt units = = =  k c 1). In the Feynman path integral formulation, the partition function is represented as a functional integral of the exponential of an effective action integrated over all fields; Here, C(β) is a normalization and the integration is performed over periodic boson loops and anti-periodic fermion loops. In quantum chromodynamics the fields f are the fermion, gluon, and ghost fields [7]. When including up to 2-loop corrections, the grand thermodynamic potential takes the form: where the sum is over quark flavors, (for our purposes i=up, down, and strange). Here, q 0 and g 0 denote the 0 th -order bag model thermodynamic potentials for quarks and gluons, respectively, while q 2 and g 2 denote the 2-loop corrections. In the last term BV is the QCD vacuum energy with B the bag constant. In most calculations sufficient accuracy is obtained by using fixed current-algebra masses. For this work we utilize a strange quark mass of = -+ m 95 from the Particle Data Group [12]. We also adopt a bag constant B 1/4 =165-240 MeV as this range should lead [2] to a critical temperature consistent with the range deduced in the low chemical potential limit from from lattice gauge theory [13].
With the parameters thus defined, the quark contribution to the thermodynamic potential is then given in terms of a sum of the ideal gas contribution plus a two loop correction from phase-space integrals over the Feynman amplitudes [7,10]: where the sum on i is over quark flavors, N c =3 is the number of colors, N g =8 is the number of gluons, and ( ) = + E p m p i 2 is the relativistic energy. The  N i denote the quark and anti-quark Fermi-Dirac distributions: The one-and two-loop gluon and ghost contributions to the thermodynamic potentials can be evaluated in a similar fashion to that of the quarks. where α s is the strong coupling constant. Here, one can immediately see from the pre-factors that the 2-loop contribution is comparable to the 0-order contribution. It is only suppressed by a factor of ∼2α s /π. For the application here we adopt the low-energy value of α s =0.33 [12].
For the massive strange quark the ideal gas contribution (equation (5)) can be easily integrated as described below. However, the two-loop correction (equation (4)) cannot be integrated numerically due to inherent divergences. It is, therefore, common to ignore the 2-loop correction or to approximate the two loop strange quark contribution with the zero mass limit [10].
This may, however, over estimate the strange-quark contribution as it ignores the Boltzmann factor suppression of the thermodynamic potential for quarks with finite mass. When the quark mass is relatively small compared to its chemical potential, this may be a reasonable approximation. However, this is not necessarily the case as one approaches the QCD phase transition at moderate values of the chemical potential.
In this paper, therefore, we attempt to extend this zero-mass quark approximation to finite masses and chemical potential followed by an analytic extension to the massless limit. We propose that this provides a more realistic input to the equation of state than the zero-quark mass limit for the 2-loop diagrams that is often employed. As an illustration we compute the bag-model QCD phase diagram for a finite-mass s-quark and contrast this with the massless limit.

Feynman diagrams
The Feynman diagrams for the fermions that make the 2-loop contribution to the thermodynamic potential are shown in figure 1. However, it is a well known problem that for finite mass all of the Feynman diagrams diverge as Λ 4 as the ultraviolet cutoff increases. Here, however, one can slightly circumvent this problem by introducing periodic regularization as in [7]. That is, the quark masses are individually treated with periodic cut-offs.
The 2-loop gluon Feynman diagrams that contribute to the potential due to the ghost loops are shown in figure 2.

Method of calculation
It is generally expected that for sufficiently high densities and/or temperature, a transition from hadronic matter to quark-gluon plasma (QGP) can occur [10]. Progress in lattice gauge theory (LGT) has shed light on the transition to a QGP in the low baryon chemical potential, high-temperature limit [13]. It is now believed that at high temperature and low density a deconfinement and chiral symmetry restoration occur simultaneously at the crossover boundary. In particular, at low density and high temperature, it has been found that the order  parameters for deconfinement and chiral symmetry restoration changes abruptly for critical temperatures of T c =145-170 MeV [13][14][15]. However, neither order parameter exhibits the characteristic change expected from a 1st order phase transition. Indeed, an analysis of many thermodynamic observables confirms that the transition from a hadron phase to a high temperature QGP is a smooth crossover [16,17]. However, it is expected that for higher chemical potential μ300 MeV, a critical point appears at which first order chiral transition can occur [13]. We note that the bag-model described here does not distinguish between the chiral and deconfinement transition. Hence, we do not describe the critical point deduced from LGT, but take that as a given. Moreover, the order of the transition requires a determination of the surface tension for nucleated bubbles of QGP [2,18]. However, there is considerable uncertainty in the surface tension. Hence, we assume an abrupt transition between hadronic and quark matter as in [18].
In spite of these limitations, it is nevertheless worthwhile for astrophysical applications to evaluate the impacts on the QCD phase diagram from both the zero-order bag model plus the massless two-loop contributions, along with effects of the finite s-quark mass on the 2-loop contribution. That is our goal in this paper.

Evaluation of feynman diagrams
The evaluation of the fermion Feynman diagrams shown in figure 1 can be reduced to the following the contour integral [7], One can then evaluate the bd + n n n , p q p factor using periodic conditions to regularize the function: where, the periodicity obeys, where T i is the function to regularize. Finally, using the on-shell condition,  where, + N F and -N F are the fermion distribution functions as defined in equation (6).
In the zero-temperature limit then, where, N F is also defined in equation (6) for ghosts. In the finite temperature m→0 limit we then have: which trivially reduces to equation (10).
Finally, for moderate temperature and μm we have the result of interest for this paper.
This expression obviously diverges exponentially for μ>m and diverges with temperature as T 2 . Nevertheless, it is useful for low to moderate chemical potential and temperatures up to and beyond ∼m s . Moreover, in the next section we describe an analytic interpolation between the massive and massless regimes that we propose is a more accurate representation of the true EoS than to assume massless quarks for the two-loop correction.

Results
In figure 3 we compare calculations of the magnitude of the the s-quark 2-loop contribution to the pressure as a function of chemical potential, where Here, we take a realistic s-quark mass of m s =95 MeV and select a typical core-collapse supernova temperature of T=10 MeV. Lines are drawn for calculations in the m s =0 limit (green line) and in the finite-mass periodicregularization of equation (21) (blue line). For finite s-quark mass, there is an obvious divergence at large μ/T due to the exp[2μ/T] term in equation (21). This limits the range of validity of this periodic regularization approach. Moreover, as can be seen in the example given in figure 3, for low chemical potential, μ<m s , the massless approximation over-estimates the pressure contribution due to the s-quark by more than an order of magnitude, and continues to overestimate the pressure by at least a factor of two up to μ∼1.5m s . Thus, one can only apply equation (21) in the low chemical-potential limit, i.e. μm s . On the other hand, in the large chemical potential limit, μ2m s one expects the thermodynamic potential to approach the massless limit of equation (20). So, we are faced with the dilemma that the function is well posed for μ<m s and μ2m s , but is not well defined between. Hence, we have considered various interpolation schemes.

Interpolation schemes
To interpolate between the small and large chemical potential limits of the EoS a Padé rational polynomial interpolation might seem physically well motivated. Indeed, there is a vast literature dealing with the Padé approximants and there exist many examples in physics of quantities that can be deduced by Padé approximants [19][20][21]. For example, the D-Padé method provides a means to extrapolate from the weak coupling to strong coupling limits of Hamiltonian Lattice QCD in the t-expansion [22,23]. Our desire here is to reconstruct the infinite series that connects between regimes of low to high values of μ/T. Hence, we seek guidance in the construction of that infinite series by Padé approximants.
The Padé approximation consists of a rational polynomial to produce an infinite series that is often a better approximation to a function than truncating its Taylor series, and may still work even in cases where the Taylor series does not converge. In the application of interest here, the Padé series involves well defined approximate functions in the low and high limits of the variable x=μ/T.
There is a well established technique for numerical Padé interpolation between the ranges of validity [21]. This is based upon the algorithm of Burlisch-Stoer [24] that generates rational functions through a recursion relation based upon Neville's algorithm that makes use of data that exist at the two limits of the interpolation [21]. Figure 4 illustrates the numerical Padé interpolation [21] based upon this algorithm in the range from from μ=75 MeV to μ=190 MeV. As can be seen by the red line in figure 4, the Padé interpolation unfortunately leads to a number of spurious poles in the Padé function. Indeed, spurious poles in the Padé function are a common feature of rational polynomial interpolation [25]. Although methods have been discussed to minimize them [19,20], we have found these to be difficult to implement and not particularly useful. In the present application these spurious poles can be traced to the fact that the application here differs from a simple truncated Taylor series. We have explicit functions in both the low chemical potential  m m (Equation (21)) and large  m m limit (equation (20)). Although the lower limit limit can be represented by a Taylor-Maclaurin series around x=(μ/T)=0, the upper limit is a simple quartic function. Achieving this limit requires a cancellation of many terms in the infinite Padé series in the limit of large μ/T.
Hence, it is difficult to impose a rational polynomial that asymptotes to equation (20) without introducing spurious poles in the recursion relation. Nevertheless, the correct interpolation is evident in the trend represented in the Padé approximant. The next obvious choice is that of a simple polynomial interpolation scheme. However, because of the rather steep transition between the behavior at low and high chemical potential a rather high order polynomial is needed.
The dashed line on figures 3 and 4 shows a 5th-order polynomial interpolation between the same limits as the Padé series. One can see that a polynomial interpolation based upon the same Neville recursion relation [21] follows the trend of the Padé interpolation without the introduction of spurious poles. Hence, a simple 5th order polynomial interpolation between the regimes obtains a result that is equivalent to the Padé interpolation but avoids the spurious poles.
As one more possibility, we have considered the application of a simple cubic spline between the small μ and large μ regimes [21]. This interpolation is shown on figure 5 and compared with the polynomial interpolations of figure 3. Here we see, that although a spline interpolation is easier to implement, it does not make use of the Although a 5th order polynomial interpolation is a possible means to generate the 2-loop contribution, we have found that this choice is difficult to implement in practical calculations. Indeed, when applied over many different temperatures a 5th order polynomial can lead to spurious fluctuations as serious as the poles in the numerical Padé method. Moreover, as can be seen in figures 4 and 5, there is inherent uncertainty of about ±30% due to the ambiguity in the choice of when to begin and end the interpolation between the massive and massless limit. For this reason (and for the fastest practical numerical applications) we suggest a simple logarithmic interpolation between μ=m s and μ=2m s . This gives a comparable uncertainty to that of the best polynomial fits and avoids completely the dangers of spurious poles and wiggles in the interpolation scheme. This is illustrated by a straight red line on figure 5.

Impact of present work on the QCD phase diagram
There are three aspects of the present work of physical relevance: 1) On the one hand, we have made a study of the importance of two-loop corrections to the QGP EoS. These are shown to be significant even in the massless limit; 2) On the other hand, we have discussed the fact that including the mass of the s-quark in the two-loop correction is divergent, but possible for low to moderate chemical potential by introducing periodic regularization of the relevant Feynman diagrams; and 3) We have shown, however, that there is a natural way to extrapolate from the low to high chemical potential regime that allows for the inclusion of the finite s-quark mass. Here we discuss these aspects and how they affect the QCD phase diagram, i.e. a plot of the critical temperature T c versus the baryon chemical potential for the transition between the hadron phase and quarkgluon plasma.
The construction of the phase diagram requires a model for both the QGP phase described here and the confined hadron phase. The equation of state for the hadron phase at very high chemical potential and density relevant to cold neutron stars, however, requires a detailed treatment of the nuclear-matter equation of state at high density. We will address this in a future work following the approach outlined in [26]. However, our main interest here is the effects of the s-quark mass at low to moderate chemical potential. For this purpose one can treat the hadron phase as a non-interacting gas of baryons and mesons that obey the usual Fermi-Dirac or Bose-Einstein statistics. , and the sum is over all baryon resonances listed in the particle data book [12], with g i the usual spin factor. Similarly, for mesons one can write: From the above thermodynamic potentials one can deduce the critical temperature for the sharp transition between the QCD and hadronic phases via a Maxwell construction [2]. Figure 6 shows an example of the impact of the formulation presented here on the QCD phase diagram. For this illustration, we chose a bag constant (B 1/4 =165 MeV). This leads to a high-temperature phase transition on the low end of the range (145 MeVT c 170 MeV) suggested by lattice QCD. We choose this value because the phase transition for lower values of B are more likely to be manifest in astrophysical environments. We note, that since we do not include hadronic interactions in the nuclear equation of state, and the modified bessel function expansion is only valid up to μ900 MeV, we do not compute the phase diagram beyond 900 MeV. This is adequate for our purpose.
The blue dot-dashed line in figure 6 shows the effect of the usually adopted assumption of only the zeroorder MIT Bag model [3][4][5][6]. The importance of adding the two-loop contribution in the massless limit is shown by the red dashed line. Here it is evident that the 2-loop contribution raises the critical temperature and/or chemical potential of the QCD phase transition. The reason for this is easy to understand. The 2-loop contribution to the thermodynamic potential enters with the opposite sign. The effect of the 2-loop diagrams is, therefore, to decrease the pressure in the QGP phase relative to that of the 0-order contribution. Since the pressure is lower one must go to higher temperature or chemical potential before the pressure of the QGP phase exceeds that of the hadron phase. Hence, a higher critical temperature results. Indeed, for this value of the bag constant, the 2-loop correction is required to obtain a critical temperature T c 145 MeV as suggested by LGT.

s-quark mass and the high-temperature phase transition
The black line on figure 6 shows the impact of incorporating the finite s-quark mass (95 MeV) into the 2-loop contribution via the procedure of periodic regularization and interpolation as discussed in the preceding sections. As expected, this correction mainly affects the regime of low chemical potential (μ300 MeV) for which the periodic renormalization and extension described above should be valid. Here one can see that adding a finite s-quark mass slightly decreases critical temperature. This is easy to understand as the effect of the finite mass is to decrease the 2-loop contribution from the s-quark until the chemical potential exceeds the s-quark mass. Thus, the diminished magnitude of the 2-loop contribution means that the pressure of the QGP phase can exceed that of the hadrons for a slightly lower critical temperature.
The portion of the QCD phase diagram in the regime of low chemical potential affected by the s-quark mass corresponds to the high temperatures of the early universe or perhaps in heavy ion collisions. There have been a number of studies of the QCD phase diagram and its effect on neutron stars in the context of chiral perturbation theory (e.g. [8,9]). It is expected that the finite-density QCD phase transition is a first-order chiral transition that occurs for baryon chemical potential from μ∼300-900 MeV. During a core collapse supernova the central core approaches and exceeds nuclear matter density with μ∼900 MeV, while the core temperature can be in excess of 50 MeV. As the core experiences a first order phase transition the baryon pressure support diminishes (although some electron pressure support remains). The initial approach to nuclear matter density induces a shock as the EoS stiffens [27]. However, the subsequent onset of a first-order transition leads to a second shock as the core re-collapses through the phase transition to form pure quark-gluon plasma [3][4][5][6].
During the initial collapse, the temperature is high enough that the massless limit is appropriate for the u and d quarks. The production of s quarks will be delayed by the time scale for weak decays, and will be further suppressed by adopting the EoS discussed here due to the Boltzmann factor ( ) -m T exp . This will cause the initial evolution through the first and second shock to proceed as though only two massless quark flavors are present. It is only during the subsequent cooling phase where the higher chemical potentials and production of massive s-quarks will be manifest.
The presence of negatively charged s quarks will neutralize the matter and diminish the isospin asymmetry. A large abundance of s quarks leads to color superconductivity in a color-flavor locked (CFL) phase in which u, d, and s quarks are paired in a symmetric and electrically neutral way. This will further decrease the pressure support of the core due to the decreased quark degeneracy pressure as u, d quarks convert to s quarks. The higher densities and temperatures in the core will lead to enhanced neutrino luminosity at late times which may be detectable. Similarly, the decreased pressure support of the core could lead to more compact neutron stars as deduced from the LIGO/VIRGO analysis of gravitational waves from the GW170817 event [28]. However, the two-loop correction introduced in this paper will suppress the s-quark content at low chemical potential. Depending upon chemical potential and temperature, this may inhibit the formation of the CFL phase and diminish the softening of the EoS relative to an EoS with massless s-quarks. This is consistent with the results of chiral perturbation theory in which the dynamical s-quark mass increases, thereby suppressing the formation of the CFL phase [9].

Comparison to other works
There have been a number of efforts to explore the phase diagram in the literature including that obtained from Lattice QCD calculations such as [29][30][31], chemical freeze-out from relativistic heavy ion collision data [32][33][34][35] and in recent NJL models which explore the impact of the phase transition on the physics neutron stars [36,37]. As noted above lattice gauge theory with a finite mass s-quark leads to a high temperature phase transition at a temperature in the range of 145 MeVT c 170 Mev [13,30], possibly favoring the lower end of the range of T c =154±9 as in [29]. A study of the free energy of 2+1 static quarks in LGT has also indicated that the chiral and deconfinement transitions occur at near the same temperature consistent with the bag-model description here [31]. Studies of the reconstructed phase diagram based upon the chemical freeze-out from relativistic heavy ion collision data indicate a critical temperature in the range 150 to180MeV [32][33][34][35]. This higher temperature is consistent with the effect of adding two-loop corrections and possibly indicates a need for a slightly higher value for the QCD bag constant. The calculation of the impact of this phase diagram on neutron stars would require a detailed model for the nuclear equation of state at high density. Nevertheless, the studies based upon the NLJ nuclear models are quite similar to figure 6 [36,37]. In [37] it was proposed that g-modes in neutron stars might be a means to detect the QCD transition. This conclusion does not change in the present work

Conclusion
We have developed a method for the incorporation of realistic quark masses into the two-loop contribution to the bag-model thermodynamic potential for quarks with finite mass and finite chemical potential. We show that the method of periodic-regularization for finite quark masses is stable for an interesting window of chemical potential up to near the s-quark mass, but diverges as the chemical potential exceeds the quark mass. Nevertheless, this is an improvement over treatments with massless quarks that can significantly overestimate the pressure contribution from quarks with finite mass at low chemical potential. A rigorous Padé interpolation between the finite mass and massless quark regimes was shown to be unstable to the generation of spurious poles in the interpolation scheme, however, we have shown that lower-order interpolation schemes gives a good representation of the Padé infinite series without the introduction of spurious poles. Nevertheless, in all the interpolation schemes there is an inherent uncertainty of order 30%. Moreover, we have found that it is impractical to implement a consistent polynomial interpolation for the wide range of temperatures and chemical potentials of relevance in practical applications. We suggest that a simple logarithmic interpolation between the s-quark mass and some multiple (say 2 times the quark mass) is within the uncertainty of the various interpolation schemes and much simpler to implement.
We have shown that implementing the two loop corrections leads to substantial changes in the QCD phase diagram. In particular, the two-loop corrections push the phase transition to higher temperatures and densities. For a fixed bag constant, this will tend to diminish observable effects from the QCD transition in supernovae and neutron stars. It will also cause a slight change in the QCD phase transition in the early universe. This could affect the relic temperatures of neutrinos (e.g. ν τ ) that decouple near the temperature of the QCD phase transition.