Equation of state of cold quark matter to O ( α 3 s ln α s )

Accurately understanding the equation of state (EOS) of high-density, zero-temperature quark matter plays an essential role in constraining the behavior of dense strongly interacting matter inside the cores of neutron stars. In this Letter, we study the weak-coupling expansion of the EOS of cold quark matter and derive the complete, gauge-invariant contributions from the long-wavelength, dynamically screened gluonic sector at next-to-next-to-next-to-leading order (N3LO) in the strong coupling constant α s . This elevates the EOS result to the O ( α 3 s ln α s ) level, leaving only one unknown constant from the unscreened sector at N3LO, and places it on par with its high-temperature counterpart from 2003. This is achieved by generalizing next-to-leading order gluon self-energies within the hard-thermal-loop limit from high temperatures and densities to zero temperature. We ﬁnd that including these screened gluonic contributions at N3LO yields a remarkably well-converged EOS, with essentially no renormalization-scale dependence. Finally, we perform a Bayesian estimation of the remaining unscreened contribution at N3LO and ﬁnd that the full EOS of cold quark matter at this order may show markedly improved convergence over the lower-order results.

Introduction.-Aproper understanding of the thermodynamics of dense strongly interacting matter is an outstanding problem in theoretical physics.This is in large part due to the infamous sign problem of lattice field theory [1][2][3][4][5], which renders nonperturbative lattice Monte Carlo techniques largely inapplicable in the region of large baryon chemical potential µ B and low temperatures T in quantum chromodynamics (QCD).
At densities above 40 times the nuclear saturation density, n 0 ≈ 0.16 fm −3 , a perturbative weak-coupling expansion in the strong coupling constant α s becomes applicable within the fundamental theory of QCD, due to asymptotic freedom.In this regime, it becomes possible to study the thermodynamics of cold (zerotemperature) quark matter (QM) directly using wellestablished thermal-field-theory tools [6][7][8][9].The equation of state (EOS) of high-density cold QM has in recent years received increasing attention as a robust high-density constraint [10] to be used when performing neutron-star EOS inference at lower densities [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27].In addition to this phenomenological application, there is great theoretical interest to study high-density cold QM, due to the rich physics arising from the dynamical screening of long-wavelength chromoelectric and chromomagnetic fields.These screening effects necessitate the development of an effective field theory of the long-wavelength gluonic modes, which goes beyond a fixed loop order in the weak-coupling expansion.
Such screening also occurs in a high-temperature quark-gluon plasma.However, at high temperatures, low-energy chromoelectric and chromomagnetic gluons can additionally receive large thermal occupation numbers ∼ 1/α 1/2 s or ∼ 1/α s , respectively [28].These large occupation numbers lead to poor perturbative convergence within the long-wavelength chromoelectric sector of QCD and necessitate a nonperturbative treatment of the chromomagnetic sector at high temperatures.Crucially, however, such a long-wavelength "Bose enhance-ment" is absent within high-density unpaired cold QM, and hence the thermodynamics of such a system always remains a perturbative problem.Furthermore, we also emphasize that this conclusion is not modified by any possible color superconducting phases that may exist at high densities, since the corrections to bulk thermodynamics from these effects are exponentially suppressed ∼ exp(−#/α 1/2 s ) [29][30][31].Since long-wavelength gluons are not Bose enhanced in cold QM, one could, in principle, extend the results for the EOS to even lower densities by tackling increasingly higher-order perturbative computationssomething which could lead to dramatic improvements within the aforementioned neutron-star EOS-inference setups.In this Letter, we compute the full EOS contributions from the screened gluonic sector at next-to-nextto-next-to-leading order (N3LO) in α s in unpaired cold QM, bringing the result to the O(α 3 s ln α s ) level and on par with its high-temperature counterpart from 2003 [32].This is achieved by deriving the next-to-leading order (NLO) screening corrections to long-wavelength gluon propagation in cold QM, generalizing the recent results in [33,34] to zero temperature.We find that our result shows very small renormalization-scale dependence for µ B where it is convergent, suggesting that the screened gluonic sector is under good perturbative control at high densities.
Structure of the weak-coupling expansion of the EOS.-In the context of cold QM with vanishing quark masses in the grand canonical ensemble, the EOS is given by the free energy density Ω = −p, where p is pressure, as a function of the quark and lepton chemical potentials.In astrophysical environments, the conditions of chargeneutral, three-flavor quark matter in equilibrium under the weak interactions (beta equilibrium) are the most relevant.These conditions reduce the EOS to a function of a single quantity µ B .Up to and including NLO, the EOS, which we henceforth identify as p(µ B ), is sensitive only arXiv:2307.08734v2[hep-ph] 13 Nov 2023 to "hard" quark and gluonic corrections from momenta K ∼ µ B , since the low-momentum quarks are Pauli blocked and the screened gluonic sector is phase-space suppressed: However, beginning at next-to-next-to-leading order (N2LO) in α s , which was computed in 1977 in the context of cold QM [35,36], the pressure becomes sensitive to long-wavelength, screened gluonic fluctuations.These field modes can be described via the hard thermal loop (HTL) effective theory [28,37] (or hard dense loop (HDL) [38][39][40]), which captures the physics of the "soft" momentum scale m E ∼ α 1/2 s µ B , corresponding to chromoelectric screening.At N2LO, because of the phase-space suppression, self-interactions between these screened gluons do not yet contribute.At the next order in α s , however, interactions involving these screened fluctuations must be included.In particular, there arise contributions from self-interactions among the fluctuating soft modes and also interactions between fluctuating soft and hard modes.The pressure correction arising from the self-interactions at N3LO has been computed within the leading-order HTL theory [41,42], while the pressure correction arising from interactions between the soft and hard modes involves corrections to the leading-order HTL theory.These corrections to highdensity HTL have been calculated within QED [43][44][45] at zero and high temperatures, and recently within QCD [33,34] at high temperatures.
In total therefore, the pressure of cold QM at N3LO can be written in the form where the three terms on the right-hand side come from the different combinations of soft and hard momentum scales.The first term p s is the purely soft contribution arising from self-interactions between screened gluon field modes.The second term p m , dubbed the mixed contribution, arises from the corrections within the HTL theory.
Finally, in addition to soft and mixed terms, there are also fully hard contributions p h , entering the expansion via four-loop vacuum graphs in full QCD (see [42] for a list).In this work, we compute the screened correction p m from the mixed sector.
The screened gluonic modes lead to nonanalyticities in the pressure, arising from logarithms of the ratio of the soft and hard scales ln(m E /µ B ).To N3LO, the pressure p of cold and dense, beta-equilibrated three-color, threeflavor (N c = N f = 3) unpaired QM with massless quarks can be cast into the following form [36,[46][47][48], (see [49] for the general N c and N f expressions) where p free = 3(µ B /3) 4 /4π 2 is the pressure of a free Fermi gas of quarks in beta equilibrium, α s = α s (Λ) is the renormalized strong coupling constant in the MS scheme at the renormalization scale Λ, and X ≡ 3Λ/(2µ B ) is a quantity which should be taken to be O(1) to minimize large logarithms.We have written logarithms of α s by grouping them in the natural expansion parameter N f α s /π, and set N f = 3.The coefficient c 3,2 of the leading logarithm was originally determined in [48].Recently, an all-order leading-logarithm resummation was conducted in [50], using this term as input.
In this work, we conclusively determine the coefficient c 3,1 (X) of the next-to-leading logarithm.It should be emphasized that this term is a well-defined and independent coefficient in the perturbative series.While the coefficient is now fully determined, the parts of the c 3,1 (X)contribution proportional to ln X could already be inferred from renormalization group invariance using lowerorder results.Likewise, given this newly determined contribution, all parts of c 3,0 (X) proportional to ln k X can now be inferred with similar arguments, thus resulting in the entire p N3LO except for one constant from the hard sector.
Overview of the calculation.-Themixed contributions correspond to the following diagrams stemming from the classification of [42] The double wavy lines correspond to soft, HTL-resummed gluons; the wavy, solid, and dotted lines correspond to hard, unresummed gluons, quarks, and ghosts, respectively; and the trace is over the suppressed Lorentz indices.Furthermore, a sum over the direction of fermionic flow is implied.In Eq. ( 3), the two shaded blobs and the respective Π's correspond to the two distinct NLO soft gluon self-energies at zero temperature: namely, the two-loop corrections given in the first pair of parentheses and the power corrections, given by the O(K 2 ) term in a small-K expansion, given in the second.The G LO is the standard one-loop resummed HTL propagator [37,38], d A = N 2 c − 1, K is the loop momentum associated with the soft gluon, and the Π's have been rescaled to be dimensionless and independent of α s , (see details and explicit expressions in [49]).In dimensional regularization and in the MS scheme, the integration measure is defined by K ≡ and Λ h is the factorization scale associated with the split between the hard and soft modes [42,51].The complete Feynman rules can be found in [42].
Note that Π 2,HTL and Π 1,Pow entering in Eq. ( 3) have been recently computed at high temperature and large µ B in general covariant gauge [33], where they have been shown to contain ln T terms that diverge in the zero-temperature limit.These terms arise because the temperature regulates certain infrared (IR) divergences associated with the hard internal gluon lines in the self-energy diagrams.The corresponding zero-temperature self-energies are computed in [49] in terms of d-dimensional integral expressions.The zerotemperature limit is achieved by using the exact integral expressions of these self-energies and isolating divergent bosonic integrals that vanish at exactly zero temperature in dimensional regularization.This procedure effectively converts, after renormalization, the ln T terms into IR 1/ε terms in the strict zero-temperature expressions.We find that the final expressions for the self-energies depend linearly on the gauge parameter ξ in a general covariant gauge, in contrast to the quadratic dependence found at high temperatures in [33].
With the renormalized self-energies in hand, the mixed contribution to the pressure can then be obtained by computing the integral in Eq. ( 3).The integral over the soft momentum K contains ultraviolet (UV) divergences that arise because the HTL theory differs from full QCD in the UV.These UV divergences must cancel with corresponding IR divergences in the hard theory, contained in p h .This cancellation has been explicitly shown in QED, where the soft sector trivially vanishes, [51], but due to the added complexity of QCD, we do not show this explicitly-however, the cancellation must occur as long as HTL does indeed correctly describe the soft physics at this order.Importantly, upon computing the radial K integral in Eq. ( 3), we find that only the specific combination of self-energies Π 2,HTL − Π 1,HTL Π 1,Pow appears.Here Π 1,HTL denotes the standard one-loop HTL self-energy [37,38], appearing in G LO .This particular combination is explicitly gauge-invariant, guaranteeing the gauge invariance of the mixed contribution to the pressure.We remark here that this should be the case since the 2-loop HTL pressure corresponding to the soft sector is known to be gauge independent [42,52], and the 4-loop hard-sector diagrams are likewise known to be algebraically gauge invariant [53].We note also that the above self-energy combination also naturally appears when one rescales the HTL effective Lagrangian including the 2-loop and power corrections to bring the kinetic term to a canonical form [34,54].
Results and discussion.-Uponcomputing the renormalized p m in dimensional regularization (details are given in [49], in general N c and N f ), we find it to have a similar form to the p s computed in [41,42], namely, where the coefficients p m i denote terms in the ε expansion and do not depend on the coupling.The factor (m E /Λ h ) −2ε arises from the integral over the softtheory loop momentum K in Eq. ( 3), while the factor (µ B /(3Λ h )) −2ε arises from the hard-theory calculation leading to the self-energies.We find the new coefficients p m i in Eq. ( 4) to be We note that this now confirms through an explicit calculation a prediction made in [42]: We then combine this result with similar expressions for p s and p h in Eq. (1) (whose expressions are given in [49]) and obtain the renormalized result for the full N3LO pressure As mentioned above, we here assume the intermediate IR and UV divergences between the different sectors to fully cancel, which in turn ensures that the Λ h dependence from the different sectors cancels.From renormalizationscale independence of the partial result (see, e.g., [51]), we are additionally able to determine the full X depen- dence of p h 0 (X), leading to the form Equations ( 4)-( 7) are our main result, and they fix the coefficients in the pressure in Eq. ( 2) to be those given in Table I, with c 0 the remaining unknown constant from the hard sector.
In Fig. 1 we show the partial N3LO pressure, neglecting only the finite hard contribution p h 0 (X) at this order.In the figure, the uncertainty of the truncation of the perturbative series is estimated by varying of X ∈ [1/2, 2], shown as a shaded uncertainty band.In this and all subsequent figures, we use the three-loop beta function when computing the running α s (Λ).We see from Fig. 1 that the pressure contribution when including the screened gluonic sector at N3LO is remarkably well converged.In particular, it has nearly vanishing renormalization-scale dependence for all µ B > 2 GeV.In fact, we have verified that one obtains a similarly well-converged result when neglecting the hard contribution at N2LO as well.This is consistent with the observation made in [41,42]: In cold QM, it is the hard terms p h (X) that drive the inevitable breakdown of perturbation theory, in stark contrast to the high-temperature case where the soft modes are responsible for the breakdown.
We turn now to an analysis of the full N3LO pressure, which we have now determined up to a single unknown constant from the hard sector c 0 .Upon substituting in different values, we find that the N3LO pressure strongly depends on c 0 .An approximate value for this constant can only be determined by computing many hard fourloop diagrams-a project which will take several years of effort.Hence, we are motivated to find an alternative way to estimate its value in order to quantify how well converged we may expect the full N3LO result to be.
Since the perturbative series for the pressure of cold QM is well-behaved up to N2LO, we turn to a Bayesian model to identify which values of c 0 are most consistent with lower-order results [55].In particular, we use the abc model of convergent series as presented in [56] and implemented in the publicly available MiHO code [57].This Bayesian model, as well as the geometric model from [58] was recently applied to the high-density perturbative-QCD EOS in [59].These models assume that the perturbative series can be approximated by independent draws from a statistical model of convergent series.Upon conditioning these models using the NLO and N2LO results FIG. 1.The N3LO pressure normalized by the free pressure as a function of baryon chemical potential, including all contributions but the hard contributions p h 0 (X) in Eq. ( 7).The shaded region shows the usual scale variation of X ∈ [1/2, 2], while the solid line is the central scale choice X = 1.using Bayes's theorem, they provide posterior distributions for the N3LO pressure as a function of µ B and X, quantifying the degree to which a given p N3LO is consistent with the lower-order results.We choose to perform the following analysis at a fixed µ B = 2.6 GeV, which is the canonical value down to which the N2LO results are typically used [11].We have checked that our conclusions remain similar if we choose a slightly different matching µ B (see also [59]).
As the N3LO pressure is now a function of a single parameter c 0 , posterior distributions of p N3LO at fixed X can be converted to distributions of c 0 .To combine the resulting c 0 distributions for different X values, we choose to marginalize over the fictitious renormalization parameter X ∈ [1/2, 2] following the procedure introduced in [58].This marginalization weights different values of the renormalization parameter in proportion to how well the abc model is converged at that X.This results in larger values of X receiving slightly larger weights, and furthermore leads to a slightly more conservative distribution than would result from the alternative scale-averaging procedure [56] in this case.
The above analysis leads to the probability distribution P (c 0 ) shown in the left panel of Fig. 2. The distribution is rather broad, with a 68% credible interval corresponding to c 0 ∈ [−29, −6], denoted in Fig. 2 by the dashed gray lines.We find that the distribution takes its highest value at c 0 = −23, denoted by the blue point in this figure.The resulting N3LO pressure corresponding to this most-consistent value of c 0 is shown in the right panel of Fig. 2, where the shaded region denotes the standard scale variation of X ∈ [1/2, 2].We find that for this value, the full N3LO pressure of cold QM lies inside the N2LO result and is well converged.In particular, the errors reach ±25% at µ B = 2.2 GeV, corresponding to about n = 27n 0 , denoted in the right panel of this figure by the dashed gray line.This is to be contracted with n ≈ 40n 0 for the N2LO result at µ B = 2.6 GeV.We note that between these densities, the relative importance of the pairing contribution to the pressure increases only by ∼ 2, so pairing effects do not grow appreciably.We thus see as an outcome of this Bayesian modeling that significant improvement to the cold-QM EOS may occur by computing the remainder of the N3LO pressure.
Conclusions.-In this Letter, we have fully computed the contributions from the screened gluonic sector at next-to-next-to-next-to leading order in the strong coupling α s .This elevates the perturbative-QCD equation of state of cold quark matter to the O(α 3 s ln α s ) level, leaving only one constant from the hard sector left to be computed, and finally places it on par with its hightemperature counterpart.We have achieved this by deriving the next-to-leading-order corrections to soft gluon propagation within a high-density medium, extending existing hard-thermal-loop results to the zero-temperature limit.
A natural follow-up to this work is the evaluation of the remaining hard constant.This would allow for the use of the full N3LO equation of state in applications, such as in studies of the neutron-star-matter equation of state.As this is an involved undertaking, it will require the development of new tools and methods.Some work in this direction has already been performed [60][61][62] and more is ongoing.

Supplemental material
Here we discuss some technical details of the computation of the gluon self-energy and the pressure computation presented in the Letter.We start by discussing the general structure of the soft gluon self-energy, to fix our notation, then present the new next-to-leading order (NLO) results for the soft gluon self-energies at large chemical potential and zero temperature T .We then continue to the calculation of the mixed diagrams, presented in the main text in Eq. ( 3).Finally, we discuss the general structure of the next-to-next-to-next-to-leading order (N3LO) pressure, including the cancellation of divergences between the different sectors, and we present results for general N c N f quarks with equal chemical potentials µ q 1 The results shown can be related to those of the main text simply by substituting µ q = µ B /3.

Gluon self-energies
At zero temperature and large density, the gluon self-energy tensor Π ab µν ≡ Π µν δ ab can be expanded in the QCD coupling α s in the schematic form where the subscript counts the order in the naive loop expansion, and Π µν resum encompasses all contributions involving (hard-thermal-loop) HTL-resummed diagrams, which we will ignore here. 2The one-and two-loop self-energies can further be expanded in small external Euclidean momentum K = (K 0 , k) as where K ≡ K/|K| and the self-energy coefficients Π µν i,q have been made α s -independent and dimensionless by introducing the in-medium effective mass parameter m E , which reads for N f identical quark flavors with equal chemical potentials µ q in d ≡ D − 1 = 3 − 2ε dimensions [42,63] with Λ h as the factorization scale and γ E as the Euler-Mascheroni constant.The gluon self-energy is a symmetric rank-two tensor spanned by four basis tensors (for a comprehensive overview see [33]).However, for the small-K coefficients shown in Eq. ( 9), only two of the four respective components turn out to be non-vanishing at T = 0, namely the three-dimensionally transverse and longitudinal components, here denoted with subscripts T and L. The respective basis tensors are given by the projectors where P µν ( K) = δ µν − Kµ Kν , i, j = 1, . . ., d, and k i being the components of k.Further, the T and L components of the small-K coefficients depend only on the Euclidean polar angle Φ through a variable x defined as x ≡ tan Φ = |k|/K 0 .With the conventions laid out above, the previously known one-loop HTL self-energy components [37,38], expanded to O(ε 2 ), read3 where we have used the shorthands and Here Li s and ζ are the standard polylogarithm and (Riemann) zeta functions respectively.We also note that the standard one-loop resummed HTL-propagator [37,38] (still with the trivial color indices suppressed) reads with ξ being the gauge parameter.
The self-energies Π µν 2,HTL and Π µν 1,Pow have recently been evaluated at nonzero T and µ q in general covariant ξ-gauge [33] (and in Feynman gauge ξ = 1 in [34]).Here, we present new results for these gluon self-energies at exactly T = 0, as those in [33] contained problematic terms proportional to ln T , which have now been converted into infrared 1/ε divergences.In order to reach the T = 0 limit, one has to consider the intermediate integral expressions leading to the self-energy: The problematic terms arise from bosonic integrals of the form ∞ 0 dp p d−4 n B (p), with n B the standard Bose-Einstein distribution function.They lead to infrared divergences proportional to (T 2 ) ε if the limit T → 0 is taken after integration, but vanish due to being scale-free when computed at T = 0.In order to get the correct T = 0 self-energies, we isolate these integrals and set them to zero, leading to a final expression that is free from ln T -type singularities.Therefore, the d-dimensional bare result for the two-loop self-energy reads Π 2,HTL L,bare = λ 1 + where we have defined the prefactor λ ≡ (e γE/2 Λ/µ q ) 2ε /Γ( d−1 2 ) and the integral which can be expressed in a closed form in terms of the Gauss hypergeometric function.The bare result becomes finite upon renormalization and can be expanded in ε yielding: For the d-dimensional bare power correction, we have: The corresponding renormalized and ε-expanded expressions read where ξ ≡ 1 − ξ, so that ξ = 0 corresponds to Feynman gauge and ξ = 1 corresponds to Landau gauge.We note that there are no terms proportional to C A present in the O(ε 0 ) and O(ε) contributions to Π 1,Pow .This occurs because at T = 0 all gluonic diagrams contributing to the power correction are scaleless and therefore vanish in dimensional regularization.The exception to this is the renormalization counterterm, which shows up as the vacuum divergence ∝ C A /ε present above.

Mixed diagrams
In this subsection, we detail the evaluation of the mixed diagrams in Eq. ( 3) using the self-energies in Eqs. ( 21), ( 22), ( 25) and (26).We start by writing Eq. ( 3) in Euclidean spherical variables as4 with the integration measure The radial integral in Eq. ( 27) can be performed by using the result for I ∈ {T, L}, yielding where we use the angle brackets ⟨•⟩ d to denote a normalised angular average in d spatial dimensions.In particular, we are interested in functions f which depend only on the angle Φ, in which case Note that the radial |K| integral resulted in the combination Π 2,HTL − Π 1,HTL Π 1,Pow in which the (modified) gauge parameter ξ cancels out, demonstrating the explicit gauge invariance of the mixed sector.The remaining angular average expanded in ε yields Tr Π 1,HTL 1−ε Π 2,HTL − Π 1,HTL Π 1,Pow 3−2ε where the leading 1/(2ε)-coefficient is given by the following analytical result The remaining two subleading coefficients contain numerical parts as well, reading and Structure of the pressure As explained in the main text and in extensive detail in [41,42], the N3LO pressure can be organized into three distinct sectors-soft, hard, and mixed-according to the kinematics of the two independent gluon lines as in Eq. ( 1).The soft sector arises from a kinematic region where both gluons are soft, diagrammatically corresponding to two-loop diagrams in the HTL theory.In dimensional regularization, with d = 3 − 2ε spatial dimensions, the result for the soft sector has the form where the coefficients p s i can be read off from [41] and for the sake of completeness they are also given below. 5On the other hand, the hard sector stems from the region where both gluons are hard, corresponding to four-loop diagrams within the naive loop expansion.The renormalized result reads schematically as with the coefficients p h i currently unknown.Finally, the mixed sector corresponds to the region where one of the gluons is soft while the other remains hard giving partially HTL-resummed three-loop diagrams (see the main text).The renormalized mixed result has the form As a result of this Letter, the coefficients p m i are given in the main text in Eq. ( 5) for N c = N f = 3 and below for general N c and N f .
With the angular average in Eq. ( 32) known to O(ε), it can be inserted into Eq.( 30), which, expanded in ε and truncating at the constant order, yields the mixed coefficients p m i .Note that in order to write Eq. ( 30) in the form

FIG. 2 .
FIG. 2. (Left) The posterior distribution of c0 within the abc model after marginalizing over X ∈ [1/2, 2] [56] at µB = 2.6 GeV (see main text).The dashed lines correspond to the 68% credible interval, and the solid point is the maximum of the distribution.(Right) The resulting N3LO pressure for most-consistent value c0 = −23.The dashed line at µB = 2.2 GeV denotes the point where the N3LO result possesses ±25% errors, and the gray lines denote contours of constant number density.

TABLE I .
List of numerical values for the coefficients appearing in Eq. (2).