Equation of state of cold and dense QCD matter in resummed perturbation theory

We discuss the Hard Dense Loop resummation at finite quark mass and evaluate the equation of state (EoS) of cold and dense QCD matter in $\beta$ equilibrium. The resummation in the quark sector has an effect of lowering the baryon number density and the EoS turns out to have much smaller uncertainty than the perturbative QCD estimate. Our numerical results favor smooth matching between the EoS from the resummed QCD calculation at high density and the extrapolated EoS from the nuclear matter density region. We also point out that the speed of sound in our EoS slightly exceeds the conformal limit.


I. INTRODUCTION
A reliable estimate of the equation of state (EoS) of cold matter at high baryon density is a vital challenge in theoretical nuclear physics. In various circumstances such as the neutron star cores, the neutron star mergers emitting gravitational waves, the supernova explosion, and the heavy-ion collisions to scan over the phase diagram of matter made out of quarks and gluons (see Ref. [1] for a review on the present status and the future direction of the heavy-ion collision), the EoS is an indispensable input for theoretical studies. Conversely, experimental data available from these extreme environments provide us with useful constraints on possible EoSs, so that some theoretical scenarios can be excluded/accepted. The most well-known and successful example along these lines is the establishment of twosolar-mass neutron stars [2], which disfavors scenarios leading to soft EoS; namely, it is unlikely for dense matter to accommodate a strong first-order phase transition [3] nor condensations of exotic degrees of freedom.
The most advanced first-principles approach from the fundamental theory of the strong interaction, i.e., quantum chromodynamics (QCD) is the lattice Monte-Carlo simulation, but the notorious sign problem ruins the importance sampling algorithm for matter at finite baryon density. Still, in parameter space where the lattice-QCD simulation is at work, the validity of alternative theoretical approaches has been tested. In particular, the Hard Thermal Loop perturbation theory (HTLpt) is the most promising resummation scheme [4][5][6][7][8] that confronts the lattice-QCD results at high temperature T . The purpose of this work is to quantify the resummation effects on the EoS of cold and dense quark matter at high baryon density n B or the energy density ε.
To sharpen novelties in our work, let us briefly summarize what has been understood so far. Since the seminal works of Refs. [9,10], we had to wait for about three decades until the perturbative QCD (pQCD) EoS was augmented with the strange quark mass M s = 0 and applied to the neutron star phenomenology [11,12], where they found that the strange mass effect is crucial. The obstacle in utilizing the pQCD EoS in neutron star physics was found to be too large scale variation uncertainty in the intermediate density region (i.e., denser than the nuclear terrain but not dense enough to justify pQCD) and the theoretical efforts are progressing toward further higher-order calculations [13,14] with hope for better convergence (see also for Refs. [15,16] for an alternative approach based on the renormalization group optimization method).
From the success of HTLpt at high T , it is a natural anticipation that the same machinery of resummation would cure the convergence problem at high baryon density or large quark chemical potential µ as well, which may reduce the scale variation uncertainty. Indeed, the parallelism between the high T and high µ cases has been established based on the transport equation approach in Ref. [17]; the high-density counterparts of HTLs are called Hard Dense Loops (HDLs). As long as a resummation prescription in the quark sector is concerned, more simply, we can just take the T → 0 limit of HTLpt to introduce "HDLpt" as considered in Ref. [18] (see also Ref. [5], and we note that the term "HDLpt" was first introduced in Ref. [19]). The HTL approximation usually neglects the bare quark mass and only the screening masses of quarks enter expressions used in Refs. [5,18]. Later on, extensive discussions about the EoS and the quark star properties have been addressed in Ref. [19]. As seen in Fig. 2 of Ref. [19], however, the HDLpt hardly remedies the convergence problem associated with uncertainty of the scaleΛ = µ − 4µ in the running coupling constant α s (Λ). In the present work, as in Ref. [12], we will employ the two-loop formula; and we will take Λ MS = 378 MeV throughout, following Ref. [12]. Previously, the absence of the bare quark mass significantly simplified technicalities as well as the realization of the β equilibrium. With equal amount of u, d, and s quarks (that is automatically the case if their masses are all neglected), the electric charge neutrality follows as it is. For quantitative descriptions of the neutron star phenomenology, however, we need to take account of the strange quark mass and solve the β equilibrium condition.
There seems to be a long way left, but the phenomeno-arXiv:2011.10891v3 [hep-ph] 14 Mar 2022 logical analyses are in need of the QCD-based EoS usable for the neutron star observables. In fact, on top of extrapolated EoSs from the nuclear side, the Bayesian analysis has been recognized as a powerful instrument for the inference analysis to identify the most likely EoS based on the observational data [20][21][22] (see Ref. [23] for a review). Recently, the Machine Learning technique has been also advocated as a complementary method to infer the EoS [24][25][26]. It would be of utmost importance to make a direct comparison of the inferred EoS candidates and the QCD-based estimates. To this end, we are urged to reduce uncertainty and widen the validity region of the pQCD or HDLpt calculations.
In this work we will report the first successful attempt to construct an EoS with smaller uncertainty from the HDLpt framework incorporating the strange quark mass effect. From the technical point of view, we adopt the resummation schemes in the gluon sector as prescribed in Ref. [4] and in the quark sector as in Ref. [18] with our own extension to cope with the strange quark mass. Our expressions are given in the form of exact integrations without any expansion in terms of the screening mass as in Ref. [7]. This paper is organized as follows: In Sec. II, we present our central results, namely the reduction of the scale dependent uncertainty in the perturbative calculation with resummation. Sec. III and Sec. IV show the calculation of the EoS within the HDLpt. In Sec. V, we show the numerical results on the speed of sound, and we take into account the order α s correction. Finally, Sec. VI summarizes this paper.

II. CENTRAL RESULTS
Since technical details are cumbersome, we shall first present our central results in Fig. 1 and then proceed to technical details later. Not to make the comparison on the figure too busy, we chose only a few representative EoSs from the nuclear side; namely, the EoS extrapolated from the chiral Effective Field Theory (χEFT) calculation [28] by the green band, the Neural Network output in the machine learning analysis [25] by the red band, and the Akmal-Pandharipande-Ravenhall (APR) EoS [29] shown by the dashed line.
The orange band in the region, ε > 10 3 MeV/fm 3 , represents the results from pQCD [12] for which we utilize the concise formula as given in Ref. [27]. Higher-order corrections could be added, but the uncertainty band is not much changed from Ref. [12]. The uncertainty band width abruptly diverges, from which it has been said that pQCD is reliable only at extreme high densities far from reality. At a glance, indeed, we should understand how difficult it is to make a robust interpolation between the nuclear and the pQCD EoSs. Now, a surprise comes from a blue narrow band that represents results from our HDLpt calculations. The uncertainty band is drastically reduced and the HDLpt EoS appears to be merged into the nuclear EoSs smoothly in the inter- other EoSs. The blue and the orange bands represent our results and the preceding results from Refs. [12,27], respectively, withΛ = µ − 4µ. The green band is from the χEFT [28]. The red band shows the EoS inferred from the Neural Networks in the machine learning analysis of the neutron star observation [25]. The dashed black line is the APR EoS extrapolated from the nuclear side [29].
mediate density region. It should be noted that the APR EoS overshoots ours, but this is due to a well-known flaw in the APR EoS, i.e., superluminal speed of sound which violates causality.
One may wonder what causes such a drastic difference on Fig. 1. We can qualitatively understand this from Fig. 2 (left) in which the baryon number density n B as a function of the quark chemical potential µ is plotted. Because the HDLpt sums the quark loops up, n B is the most sensitive quantity affected by the resummation in the quark sector. It is an interesting and reasonable observation that n B is suppressed at fixed µ after the resummation: thermodynamic quantities are dominated by quark quasi-particles, and in HDLpt, quark excitations are more screened by self-energy insertions, as compared to pQCD treatments. Therefore, on Fig. 1, the corresponding µ for a given ε becomes larger, and the corresponding running coupling α s (Λ = ξµ), where ξ = 1, 2, 4, is smaller. This qualitative argument partially accounts for the reduction of the uncertainty band, but not fully yet. As shown in Fig. 2, if we plot the pressure P , the baryon number density n B , and the energy density ε as functions of µ, respectively, the uncertainty bands are wider than Fig. 1. In Fig. 2 (left), we overlay a horizontal line at n B = 10 n 0 to find the values of corresponding µ for differentΛ. The values of P at these µ's are shown in Fig. 2 (right) with the same markers. Importantly, the marker for P (Λ = µ) is out of the plot range. Owing to the suppression in n B leads to the situation that P (ε) withΛ = µ and that withΛ = 4µ happen to stay close, which narrows the uncertainty band on Fig. 1. There might be a deep reason (e.g., scaling properties) for this behavior, and further investigations are in progress.
For the astrophysical application, we need P (ε) or P (n B ) rather than P (µ). The condition that P (ε;Λ) is insensitive to the scaleΛ is dP (ε;Λ)/dΛ = 0, i.e., where µ B = 3µ is the baryochemical potential. Substituting the thermodynamic relation ε = −P + µ B n B this relation reduces to In the conventional argument, the reduction of the first terms in Eqs. (1) and (2) has been the central issue, but we point out that ∂P/∂Λ = 0 is only a sufficient condition for Eqs. (1) and (2). Albeit ∂P/∂Λ = 0, the inclusion of the latter term can cancel the scale-dependence; Fig. 2 is the concrete realization of such cancellation.

III. FORMULATION
Let us explain the formulae and procedures to obtain our results in Fig. 1. Dense matter in the neutron star reaches the β equilibrium; d u + e − +ν e and s u + e − +ν e indicating the relations between quark chemical potentials as µ u = µ + 2 3 µ Q and µ d = µ s = µ − 1 3 µ Q where µ Q is the electric chemical potential. Since electrons are negatively charged, µ e = −µ Q , and we can fix µ Q from the charge neutrality, i.e., n Q − n e = 0 with n Q = ∂P/∂µ Q and n e = µ 3 e /(3π 2 ) neglecting the electron mass.
Since the most crucial extension in this work is the inclusion of the bare quark mass, we will write down the explicit expressions in the quark sector. In our notation for flavor-f quarks the bare mass is M f and the screening mass is m qf . The bare mass should be scale dependent .
In the T → 0 limit the HDLpt pressure, P HDLpt , is given by the gluon loop and the quark loop with the selfenergy insertions; namely, where ∆P g and ∆P q subtract the ultraviolet divergences. The gluon part with an appropriate subtraction by ∆P g ∝ 1/ (where the spatial dimensions are d = 3 − 2 in the dimensional regularization) is A constant, C g , is an integral over a function involving the gluon self-energy and numerically estimated as C g ≈ 1.17201 in the dimensional regularization. Here, m D is the gluon screening mass induced by µ, i.e., m 2 D ≡ (2α s /π) f µ 2 f . We note that the bare quark masses in the hard loops are neglected commonly in the HTL approximation (see Ref. [30] for a standard textbook). The gluon sector is intact, so we just refer to Refs. [4,6,7] for further details.
The quark part appears from the flavor-f quark loop, i.e., P q,f = tr ln G −1 k − M f − Σ(k 0 , k) and k 0 = iω n + µ f for flavor-f quarks withω n being the fermionic Matsubara frequency. For the selfenergy expression, Σ, we need to introduce the following notations according to Refs. [7,18] , and the flavor-f quark screening mass is m 2 The fermionic HTLpt function in d = 3 − 2 spatial dimensions is: Then, the self-energy for flavor-f quarks is expressed as In this work, we neglect the bare quark mass dependence in Σ(k 0 , k); this treatment can be justified under the HDL approximation. In principle this effect can be taken into account by using the effective action presented, e.g., in Ref. [31]. The expression will, however, be extremely complicated, so we will simply neglect here. Nevertheless, it is unlikely that the bare quark mass plays an important role for our main results, i.e., the reduction of the scale dependent uncertainty.
The paramount advance in this work is the inclusion of bare mass M f , and the quark pressure deviates from Refs. [7,18]. Let us first write down our final expression and then explain the notations next. In the flavor-f quark sector the pressure contribution reads: We introduced C q and D q as functions of η f ≡ 1 + M 2 f /(2m 2 qf ). These definitions involve the following functions: whereω is a dimensionless and continuous variable. Then, C q and D q are given by We note that D q (η f → 1) → 0 and C q (η f → 1) ≈ −0.03653 as is consistent with Ref. [7]. The next term, P qp,f , in Eq. (7) is the quasi-particle contribution given by (12) We note that the ideal term ∝ µ 4 f is subtracted in the above expression since we doubly pick up two pole contributions at ω f ± . In Ref. [18] the quasi-particle contribution was defined by taking the m 2 qf derivative/integration, so that only the difference from the ideal term was considered by construction, and the ideal term was not subtracted but added. Here, the quasi-particle poles, ω f ± , are solutions of the following implicit equations, i.e., and Q 1 (x) ≡ xQ 0 (x) − 1 being the Legendre functions. Finally, the last term in Eq. (7) represents the contribution from the Landau damping, which reads: The integrand is given by tan θ qf = Y/X where In this case k ≥ ω holds and the argument of Q 0 should be k/ω, not ω/k. We also note that the subtraction at finite M f is mass dependent, i.e., ∆P q = m 4 qf D q (η f )/(2 ). For numerical calculations, we took M u = M d = 0 and M s (2GeV) = 100 MeV. For N c and N c in α(Λ) and M s (Λ) we took N c = N f = 3. This completes the explanation of the formulation necessary to draw Fig. 1.
In Fig. 3, we show the EoSs calculated based on the formulation presented above. In Fig. 3 (Left), we show the EoS in the form of P (ε). This is the same plot as Fig. 1 above, but with an extended region of the energy density. Because of the uncertainty out of control at lower energy density it is reasonable to truncate the plot around ε 500 MeV/fm 3 .
In Fig. 3 (Right), we show the EoS in the form of P (µ). It is evident that the scale variation uncertainty in HDLpt is not small as compared with the pQCD results. Therefore, it is a quite nontrivial discovery that the scale variation uncertainty in P (n B ) is significantly smaller than that in P (µ).

IV. DETAILS OF INTEGRATION: THE QUARK CONTRIBUTION TO THE PRESSURE
Here, we will elaborate the details of integration that appears in the derivation of Eq. (7) in the previous section. The quark part of the pressure appears from the flavor-f quark loop: where we write the sum-integral as {K} = T ωn k in d = 3 − 2 spatial dimensions for the momentum integration. The functions A 0 and A S are defined above. We note that P q,f in Eq. (17) can be regarded as a leading contribution in the 2PI or the Cornwall-Jackiw-Tomboulis (CJT) formalism [32,33]. This explains why Eq. (17) misses an additional term, tr ΣG f , that may be responsible for the deviation of O(α s ), which will be studied below. We recast the Matsubara sum into the contour integral along C as depicted in the left panel of Fig. 4. We can deform the contour C into C qp ∪ C Ld , see the right panel of Fig. 4. We identify the terms from C qp and C Ld with the quasiparticle contribution and the Landau damping contribution, respectively, according to Refs. [5,7]: The quasiparticle contribution to the integral (see the right panel of Fig. 4) is where we defined as Disc , and dropped an irrelevant infinity from the upper bound of the ω-integration. The dispersion relation for quarks ω f ± is obtained by solving Eq. (13) above. For the moment we can drop the third term in Eq. (20) that is independent of T and µ f (which will be reassembled later). Finally, we obtain: which completes the derivation of Eq. (12) above. The s = −1 term in the sum of Eq. (20) vanishes at T → 0 because of the step function θ(−µ f − ω f χ ). The Landau damping contribution to the integral is In the last line we introduced [with X and Y defined in Eqs. (15) and (16) above, respectively]: Again, we only keep the T and µ f dependent parts in Eq. (22), so that the T → 0 limit leads to which completes the derivation of Eq. (14) in the previous section.
We here reassemble the T and µ f independent terms that we dropped above. To this end it is convenient to think of the T = µ f = 0 limit in Eq. (18), in which the Matsubara sum reduces to T n → ∞ −∞ dω 2π , so that the pressure reads: where we used the following integral: The function f ± (ω, η f ) with η f ≡ 1 + M 2 f /(2m 2 qf ) is defined as in Eqs. (8) and (9) above. The limit of → 0 gives: The constants C q and D q are defined in Eqs. (10) and (11) in the previous section, respectively. The ultraviolet divergence is subtracted by the term ∆P q in Eq. (4) above: In this way the above procedures complete the derivation of Eq. (7) in the previous section.

V. DISCUSSIONS
Here, we discuss the speed of sound that could exceed the conformal limit, and the robustness against the O(α s ) corrections to match the conventional pQCD calculation.

A. Speed of sound
The EoS from our resummed perturbation theory has a notable feature in addition to the smaller uncertainty. We have calculated the speed of sound, c 2 s = ∂P/∂ε, which is depicted in Fig. 5. To make clear the relevance to the neutron star environment, we chose the horizontal axis as the baryon number density n B in the unit of the normal nuclear density n 0 .
There is an empirical conjecture to claim that the speed of sound may not exceed the conformal limit, i.e., c 2 s = 1/3. In the high density limit, asymptotically, all mass scales and interactions are negligible and the conformal limit should be eventually saturated. In the pQCD calculation, the first correction from the conformal limit is negative, so that the conformal limit is approached from c 2 s < 1/3 with increasing density. Also at finite tem-perature, the lattice-QCD results demonstrate that the conformal bound c 2 s < 1/3 holds [34]. Known examples of QCD calculations seem to respect the conformal limit (see Ref. [35] for an exception at finite isospin chemical potential). However, no field-theoretical proof exists to guarantee c 2 s < 1/3. The recent analysis based on neutron star data, especially the two-solar-mass condition, indeed suggest a possibility of c 2 s > 1/3 at sufficiently high baryon density [36][37][38]. Figure 5 shows that our resummed EoS slightly violates the conformal bound and c 2 s approaches 1/3 from above. It is evident that our result is a counterexample to the conjecture of c 2 s < 1/3. The quantitative difference is numerically small between EoSs from our HDLpt and pQCD, and the violation of the conformal bound is tiny, but this comparison on Fig. 5 implies that one should be careful about the robustness of the speed of sound bound (see, for example, discussions in Refs. [39]).

B. O(αs) correction
The HDLpt has a deviation of O(α s ) in the pressure, as mentioned in the beginning, from the conventional pQCD  calculation. Our HDLpt predicts c 2 s > 1/3 even if we add a correction to match the O(α s ) terms. For analytical simplicity we will show the calculation in the massless case only. It is known that the expansion of P HDLpt in powers of m qf /µ f 1 gives, for N c = 3 [18]: where the ideal pressure is P ideal = N c N f µ 4 f /(12π 2 ). The conventional pQCD result is [9,10] Therefore we can match the O(α s ) terms by adding the following correction to P HDLpt : In Fig. 6 we plot the speed of sound evaluated by P HDLpt and P HDLpt + P corr both in the massless case. The figure 6 clearly shows that even with the P corr correction, the speed of sound still approaches c 2 s = 1/3 from the above as the density increases. This implies that c 2 s > 1/3 could be attributed to the higher order effects from the resummation.

VI. SUMMARY
In this work we showed results with the smaller scale variation uncertainty for the cold dense matter EoS in the form of P (ε). The formalism we adopted here is the HDLpt, which has already been successful in finite temperature QCD. The important observation is that, as compared to the pQCD calculation, quarks are screened by self-energy insertions, and the baryon density is suppressed. This means that the corresponding chemical potential for a given baryon density is shifted to be larger, particularly forΛ = µ. It was the source of the large uncertainty in the pQCD calculation, so the improvements forΛ = µ helps lessen the uncertainty band. We also emphasize the importance of the inclusion of a bare quark mass and we numerically solved the β equilibrium and charge neutrality conditions. Our treatments with the bare quark mass are messy, but contributions from finite strange quark mass are crucial for the realistic environments of neutron stars under the β equilibrium. Our results constitute a QCD-based example of the conformal limit violation at finite density, which can be in consonance with the state-of-the-art neutron star observations. It would be an exciting program to apply our EoS to the neutron star phenomenology. We will report phenomenological implications soon.