# All order resummed leading and next-to-leading soft modes of dense QCD pressure

Loïc Fernandez, Jean-Loïc Kneur

► To cite this version:

Loïc Fernandez, Jean-Loïc Kneur. All order resummed leading and next-to-leading soft modes of dense QCD pressure. Physical Review Letters, 2022, 129 (21), pp.212001. 10.1103/Phys-RevLett.129.212001. hal-03347966

HAL Id: hal-03347966

https://hal.science/hal-03347966

Submitted on 13 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# All Order Resummed Leading and Next-to-Leading Soft Modes of Dense QCD Pressure

Loïc Fernandez and Jean-Loïc Kneur

Laboratoire Charles Coulomb (L2C), UMR 5221 CNRS-Université Montpellier, 34095 Montpellier, France

<figure><img src="image_2.png" /></figure>

(Received 22 September 2021; revised 9 June 2022; accepted 21 October 2022; published 17 November 2022)

The cold and dense QCD equation of state at high baryon chemical potential $\mu_B$ involves at order $\alpha_s^2$ an all-loop summation of the soft mode $m_E \sim \alpha_s^{1/2} \mu_B$ contributions. Recently, the complete soft contributions at order $\alpha_s^3$ were calculated using the hard thermal loop formalism. By identifying massive renormalization group properties of the hard thermal loop theory, we resum to all orders $\alpha_s^p$, $p \ge 3$ the leading and next-to-leading logarithmic soft contributions. We obtain compact analytical expressions that show visible deviations from the state-of-the art results, and noticeably reduce residual scale dependence. Our results should help to reduce uncertainties in extending the equation of state in the intermediate $\mu_B$ regime, relevant in particular for the phenomenology of neutron stars.

DOI: 10.1103/PhysRevLett.129.212001

**Introduction.**—At sufficiently high temperature or density, the asymptotic freedom property of quantum chromodynamics (QCD) naively provides a weak coupling perturbation theory (PT) to address the quark-gluon plasma physics. However, severe infrared (IR) divergences spoil a naive PT approach, giving poorly convergent results at successive orders, unless at extremely high temperatures and/or densities (see, e.g., [1] for reviews). In contrast, today's powerful lattice simulations offer an alternative numerical nonperturbative *ab initio* solution. So far lattice simulations have been very successful in the description of the QCD crossover transition at finite temperatures and near vanishing baryonic densities, with results [2] also relevant to confront the experimental data from heavy ion collisions in this specific region of the phase diagram. However, the notorious sign problem [3] prevents simulations at high densities, equivalently high chemical potential $\mu_B$, to explore the more complete QCD phase diagram and in particular at $\mu_B$ values pertinent to the physics of neutron stars [4]. Alternatively, more analytical approaches to resum thermal and in-medium PT have been developed and refined over the years (see, e.g., [1,5–9]), improving the bad convergence generically observed even at moderate coupling values. More recently, a renormalization group (RG) resummation approach giving sizeably improved renormalization scale uncertainties was developed [10,11].

In the following, we focus on cold and dense QCD, $T = 0$, $\mu_B \neq 0$, which implies a number of simplifications. For strongly coupled matter at high baryonic density, the

dynamical screening of color charges manifests in a screening mass defining a soft scale $m_E \sim \sqrt{\alpha_s \mu_B} \ll \mu_B$, where $\alpha_s$ is the QCD coupling. This is linked to IR divergences in the naive perturbative pressure, to be appropriately resummed, resulting in nonanalytic $\ln \alpha_s$ dependences in the perturbative expansion. This phenomenon occurs first at order $\alpha_s^2$, giving a contribution $\alpha_s^2 \ln \alpha_s$ established for massless quarks long ago [12]. The pressure was extended much later to massive quarks [13–15], but for a very long time no higher order was available. These soft terms are not the full contributions at $\alpha_s^{p \ge 2}$ orders, and should be completed by hard contributions $\mu^4 \alpha_s^p$ calculable from standard PT, known exactly to date up to order $\alpha_s^2$ [12,16]. Yet, the soft terms constitute a well-defined subset, relevant for the convergence of the weak coupling expansion. The combined hard and soft contributions exhibit sizeable residual dependence in the (arbitrary) renormalization scale (although less severe than for thermal QCD), leading to systematic uncertainties. It is therefore crucial to push further the weak coupling expansion, as one expects that higher perturbative orders may reduce uncertainties in the equation of state (EOS), in particular in a regime presumably relevant to the physics of neutron stars [4]. Given the abundant new data on compact stellar objects from astrophysics, and the rapidly developing interplay between gravitational wave physics, QCD, and nuclear calculations, it is timely to try to further tighten the existing gap in the moderate $\mu_B$ regime between the reliable perturbative QCD at high $\mu_B$ and reliable EOS in the low $\mu_B$ nuclear regime. Recently, the complete soft terms at the next order, $\alpha_s^3 \ln^p(m_E)$, $0 \le p \le 2$, were obtained in [17–19] from involved calculations using the hard thermal loop (HTL) formalism to unprecedented $m_E^4 \alpha_s$ order.

Our main purpose in this Letter is to calculate and resum the soft leading logarithms (LLs) and next-to-leading

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Funded by SCOAP³.
logarithms (NLLs) to all $\alpha_s^p$ orders, $p \ge 3$. As particular case it provides an independent simpler derivation of the first LL term $\alpha_s^3 \ln^2 \alpha_s$, consistent with [19].

**Cold quark matter state-of-the-art pressure.**—For $N_f$ massless quarks with $\mu_q = \mu \equiv \mu_B/3$ the presently known quark matter weak expansion pressure reads

$$P^{\text{cqm}} = P_f \left[ 1 - \frac{2}{\pi}\alpha_s - \frac{N_f}{\pi^2}\alpha_s^2 \ln \alpha_s - 0.874355\alpha_s^2 - 2d_A \frac{(11N_c - 2N_f)}{3(4\pi)^2} \ln \left(\frac{M_h}{\mu}\right) \alpha_s^2 \right] + P_{\alpha_s^3}^{\text{soft}}, \quad (1)$$

with $d_A = N_c^2 - 1$ ($N_c = 3$), and other terms specified and commented on below. The leading order (LO) massless quark loop gives the free gas pressure $P_f = N_c N_f \mu^4 / (12\pi^2)$, and next-to-leading-order (NLO) gluon exchange gives the $O(\alpha_s)$ term. At $\alpha_s^2$ order the three-loop graphs in standard perturbation give [12,16] the hard contributions $\sim \alpha_s^2 \mu^4 / P_f$ in Eq. (1), also involving an arbitrary renormalization scale $M_h$. The nonanalytic $\alpha_s^2 \ln \alpha_s$ term arises from resumming [12] the set of all order soft contributions, the “ring” graphs (on left in Fig. 1). A well-known important feature is that the all-loop summation gets rid of initially IR divergences, a mechanism related to the dynamical gluon screening mass. At $T = 0$, $\mu \neq 0$, $m_E$ is obtained from the sole quark-loop contribution to the self-energy:

$$m_E^2 = 2 \frac{\alpha_s}{\pi} \sum_f \mu_f^2 = 2 \frac{\alpha_s}{\pi} N_f \mu^2. \quad (2)$$

Accordingly, in a rough picture, the all-loop soft mode summation can be essentially obtained from a one-loop calculation, within an alternative framework with a massive gluon (the graph on the right in Fig. 1), giving a contribution $P \sim m_E^4 \ln m_E \sim \mu^4 \alpha_s^2 \ln \alpha_s$. This picture is rigorously embedded in the HTL formalism [7,20] that essentially provides a gauge invariant effective field theory (EFT) consistently including all HTL contributions with dressed momentum-dependent self-energies and vertices, and involving the screening gluon mass. The $\alpha_s^3$ last contribution in Eq. (1) was calculated [17,19] from HTL graphs with only gluons, NLO corrections to the right graph in Fig. 1. The $m_E^4 \alpha_s \ln^2(m_E) \sim \alpha_s^3 \ln^2 \alpha_s$ was first obtained in [17], and recently the remnant soft terms were completed, obtaining [18,19]

<figure><img src="image_3.png" /><figcaption>FIG. 1. Connection between ring sum and one-loop HTL.</figcaption></figure>

$$P_{\alpha_s^3}^{\text{soft}} = \frac{N_c d_A \alpha_s m_E^4}{(8\pi)^2} \left( \frac{p_{-2}}{4\epsilon^2} + \frac{p_{-1} - 2p_{-2} \ln \frac{m_E}{M_s}}{2\epsilon} + 2p_{-2} \ln^2 \frac{m_E}{M_s} - 2p_{-1} \ln \frac{m_E}{M_s} + p_0 \right), \quad (3)$$

$$p_{-2} = \frac{11}{6\pi}, \qquad p_{-1} \approx 1.50731(19), \qquad p_0 \approx 2.2125(9). \quad (4)$$

All the $\mu$ dependence arises solely from the screening electrostatic mass, $m_E$ from Eq. (2). Notice in Eqs. (1) and (3), besides $M_h$, the different scale $M_s$ introduced in [18] (our $M_s$ scale corresponds to $\Lambda_h$ in [18]), for the soft sector. As explained in [18,19], the presently unknown $\alpha_s^3$ hard (and mixed soft-hard) contributions are expected to cancel the remnant UV divergences and soft $\ln^p(m_E/M_s)$ in Eq. (3), to let only $\ln^p \alpha_s$ and $\ln^p(M_h/\mu)$ terms ($p = 1, 2$). In the absence of such explicit cancellations, since the soft terms can be treated as a separate $m_E$-dependent sector, to avoid large logarithms it appears sensible to choose [18] $M_h \sim O(\mu)$ and $M_s \sim O(m_E)$.

**HTL one-loop pressure.**—Our starting expression is the one-loop HTL pressure, calculable from the graph in Fig. 1. In the present cold quark matter context it is evaluated at $T = 0$, with (yet unspecified) mass $m_g$, obtaining

$$P_{\text{LO}}^{\text{HTL}} = \frac{d_A m_g^4}{(8\pi)^2} \left[ \frac{1}{2\epsilon} + C_{11} - L + \epsilon(L^2 + C_{21}L + C_{22}) \right] \\ \equiv m_g^4 \left[ -\frac{a_{1,0}}{2\epsilon} + a_{1,0}L + a_{1,1} + O(\epsilon) \right], \quad (5)$$

where $L = \ln(m_g/M)$ and $M$ is an arbitrary $\overline{\text{MS}}$ scale. In Eq. (5) we also introduced a convenient notation for later purpose. The coefficient $C_{11} \sim 1.17201$ was first obtained in [21]. As explained below, however, to work out complete NLL expressions we also need the $O(\epsilon)$ terms, not previously available to our knowledge, that we have evaluated. We obtain $C_{21} = -2C_{11}$, $C_{22} \approx 2.16753$ [22].

**RG resummation.**—We consider the RG operator (see, e.g., [23]) that parametrizes the scale variation in a (massive) theory:

$$M \frac{d}{dM} \equiv M \frac{\partial}{\partial M} + \beta(g^2) \frac{\partial}{\partial g^2} - \gamma_m^g(g^2) m_g \frac{\partial}{\partial m_g}, \quad (6)$$

where $M$ is a renormalization scale, $\beta(g^2) \equiv dg^2/d \ln M$ dictates the running coupling, and the anomalous mass dimension $\gamma_m^g(g^2) \equiv d(\ln m_g)/d \ln M$ is further discussed below. In our convention with $g^2 \equiv 4\pi\alpha_s$,

$$\beta(g^2) = -2b_0^g g^4 - 2b_1^g g^6 + \dots, \\ \gamma_m^g(g^2) = \gamma_0^g g^2 + \gamma_1^g g^4 + \dots, \quad (7)$$
where $b_0^g, b_1^g$ are the QCD “pure gauge” contributions

$$ (4\pi)^2 b_0^g = \frac{11N_c}{3}, \quad (4\pi)^4 b_1^g = \frac{34N_c^2}{3}. \quad (8) $$

Given that HTL involves a gluon mass term $m_g$, from a RG standpoint [24] we consider $m_g$ blind to its precise dynamical origin as a screening mass, motivating us to use the massive RG Eq. (6). Indeed, despite the nonlocal HTL Lagrangian, the next-to-next-to-leading-order (NNLO) perturbative HTL calculations give new $m_g$-dependent UV divergences and related counterterms having a renormalizable form, as will be seen below, thus defining EFT anomalous dimensions in standard fashion. It is important to stress that the RG coefficients for such a massive theory are $T=0$ entities by definition, even though they enter thermal or in-medium contributions as well. More precisely, within $T \neq 0$ HTL calculations, the divergences from $m_g \neq 0$ occur in two-loop order $\alpha_s(m_g^2 T^2, m_g^3 T)$ terms, and the corresponding (unique) one-loop counterterm $\Delta m_g$ was obtained first in [6]. Using this $\Delta m_g$ and the standard relation between bare mass $m_g^B$, $Z_{m_g}$ counterterm and Eq. (7): $m_g^B \equiv m_g Z_{m_g} \approx m_g[1 - g^2\gamma_0^g/(2\epsilon) + \mathcal{O}(g^4)]$, we easily identify

$$ \gamma_0^g = \frac{11N_c}{3(4\pi)^2} \equiv b_0^g. \quad (9) $$

Namely, the LO interactions from HTL that contribute to renormalize $m_g$ give a divergent contribution identical to the one defining $b_0^g$. Although striking, this equality of pure gauge $b_0^g$ and $\gamma_0^g$ is merely a one-loop order accident. Incidentally, it is worth noting that the same result [Eq. (9)] was obtained independently from a localizable, renormalizable gauge-invariant setup for a (vacuum) gluon mass [25]: this is not a coincidence since those universal RG quantities are vacuum quantities independent of $T$, $\mu$. The two-loop order $\gamma_0^g$, entering the NLO RG [Eq. (7)] in our construction, has also been calculated from the same $T = \mu = 0$ formalism, with the result [26] $(4\pi)^4\gamma_1^g = 77N_c^2/12$. Furthermore, the mass renormalization within Eq. (5), $P_{\text{LO}}^{\text{HTL}}(m_g \rightarrow m_g Z_{m_g})$, generates additional terms that combine with genuine two-loop contributions [Eq. (3)]. Importantly, the unwanted nonlocal $\ln(m_E/M_s)/\epsilon$ divergence in Eq. (3) exactly cancels in those combinations, while (local) remnant divergences after mass renormalization are renormalized by vacuum energy $\mathcal{E}_0 (\equiv -P)$ counterterms, always necessary in a massive theory. According to Weinberg's theorem [23], such local counterterms prove the renormalizability of the $T=0$ HTL pressure at NLO $\alpha_s m_g^4$, i.e., next-to-next-to-next-to-leading-order $\alpha_s^3$. The LO vacuum energy counterterm is the one determined in HTL [6,21], $\Delta\mathcal{E}_0^{(1)} = d_A m_g^4/(8\pi)^2/(2\epsilon)$. At NLO its expression is given in

[22]. Remark that renormalizing $m_g$ in Eq. (5) also modifies the “finite” coefficients in Eq. (4) as

$$ \begin{aligned} 2p_{-2} &\rightarrow p_{-2}, \\ p_{-1} &\rightarrow p_{-1} - \frac{8\pi\gamma_0^g}{N_c} \left(C_{11} - \frac{1}{4}\right) \approx p_{-1} - 0.5381, \\ p_0 &\rightarrow p_0 - \frac{8\pi\gamma_0^g}{N_c} \left(C_{22} - \frac{C_{11}}{2}\right) \approx p_0 - 0.9229. \end{aligned} \quad (10) $$

In [19] nontrivial cancellation mechanisms between soft $\alpha_s m_E^4$ and (presently not known) $\alpha_s^3$ hard contributions are convincingly argued to remove the divergent terms in Eq. (3), and to obtain the correct $\alpha_s^3 \ln^2 m_E$ coefficient, $2p_{-2} \rightarrow p_{-2}$. In our massive renormalized scheme, the right $p_{-2}$ originates directly from standard mass renormalization cancellation mechanisms. This simple alternative picture appears consistent with generic EFT constructions [27], where all EFT UV divergences are renormalized, defining EFT anomalous dimensions, sufficient to extract $\ln m_E$ dependencies from the RG.

Having identified these key RG ingredients at $T=0$, obtaining the LL and NLL at successive orders follows a well-established procedure [23]. Power counting [19] dictates the soft pressure expansion to all orders

$$ P_{\text{soft}} \sim m_g^4 \sum_{p=1}^{\infty} (g^2)^{p-1} \sum_{l=0}^{p} a_{p,l} \ln^{p-l} \left( \frac{m_g}{M} \right), \quad (11) $$

with $M$ the (MS scheme) renormalization scale. The $a_{p,0} \ln^p(m_g/M)$, $p \ge 1$ are the leading logarithmic (LL), $a_{p,1} \ln^{p-1}(m_g/M)$, $p \ge 2$ the next-to-leading logarithmic (NLL) coefficients and so on, with $a_{p,p}$ the nonlogarithmic coefficients at successive orders. Upon applying the RG Eq. (6) on Eq. (11) considering $g$ fixed, we obtain recurrence relations. First, for the LL series, we obtain

$$ -p a_{p,0} = [4\gamma_0^g + 2b_0^g(p-2)]a_{p-1,0}, \quad p \ge 2, \quad (12) $$

where $a_{1,0}$ is given in Eq. (5). Before proceeding it is worthwhile to give a first concrete outcome of Eq. (12) that immediately gives at NLO

$$ a_{2,0} = -2\gamma_0^g a_{1,0} \equiv -2b_0^g a_{1,0}, \quad (13) $$

determining the LL term $g^2 m_g^4 a_{2,0} \ln^2(m_g/M)$. Upon using Eqs. (5), (8), it gives $(4\pi)a_{2,0} = N_c d_A/(8\pi)^2 p_{-2}$, which perfectly matches Eqs. (3), (4) for $m_g \equiv m_E$ with $p_{-2}$ in Eq. (10). We stress that all LLs in Eq. (12) rely on the sole one-loop $\ln(m_g/M)$ coefficient $a_{1,0}$ in Eq. (5) and thus are completely independent from the results in [17,19], where $m_E^4 \alpha_s \ln^2 m_E$ was obtained from involved two-loop HTL calculations. To obtain the latter from the above simple RG relations the important feature is having identified $\gamma_0^g$.
More interestingly, Eq. (12) gives all higher order LL coefficients, a new result that we resum explicitly below. One obtains similarly the NLL series (defined for $p \ge 2$)

$$ (1-p)a_{p,1} = [4\gamma_0^g + 2b_0^g(p-2)]a_{p-1,1} \\ + [4\gamma_1^g + 2b_1^g(p-3)]a_{p-2,0} \\ + \gamma_0^g(p-1)a_{p-1,0}. \quad (14) $$

We recall that at a given perturbative order $g^{2p}$, the only new terms to calculate are the single logarithm $a_{p,p-1}$ and nonlogarithmic $a_{p,p}$ terms. Both the LL and NLL series above are convergent and can thus be resummed.

**Perturbatively RG invariant pressure.**—Before giving resummed LL and NLL expressions, it is convenient to develop a related important ingredient of our construction, namely to restore a RG invariant (RGI) massive pressure. Indeed, applying Eq. (6) to Eq. (5) gives a remnant scale dependence at leading HTL order, $\sim m_g^4 \ln M$, a well-known feature of massive theories. The appropriate way to deal with this is to realize that the RG Eq. (6) has an inhomogeneous term, defining a (gluon) vacuum energy anomalous dimension $\hat{\Gamma}_0^g(g^2) = \Gamma_0^g + \Gamma_1^g g^2 + \dots$, by analogy with other massive sector vacuum energies [28,29]:

$$ \frac{dP_{\text{HTL}}}{d \ln M} = -m_g^4 \hat{\Gamma}_0^g(g^2) = \frac{d}{d \ln M} \left( m_g^4 \sum_{k \ge 0} s_k^g g^{2k-2} \right), \quad (15) $$

$\hat{\Gamma}_0^g(g^2)$ being related to the vacuum energy counterterms (see Ref. [22]). Thus, an RGI combination is $P_{\text{RGI}} = P_{\text{HTL}} - m_g^4 \sum_k s_k^g g^{2k-2}$, where the $s_k^g$ coefficients are most simply determined [24] perturbatively from the second equality in Eq. (15): at LO and NLO we obtain, respectively,

$$ s_0^g = -\frac{a_{1,0}}{2(b_0^g - 2\gamma_0^g)} = -\frac{d_A}{2(8\pi)^2 b_0^g}, \quad (16) $$

$$ s_1^g = a_{1,1} + \frac{a_{2,1}}{4\gamma_0^g} + \frac{a_{1,0}}{4} + \frac{s_0^g}{2\gamma_0^g}(b_1^g - 2\gamma_1^g), \quad (17) $$

where $a_{1,1} = -a_{1,0}C_{11}$ in Eq. (5). Embedding Eq. (16) within Eq. (12) [30], one can easily resum the LL series as

$$ P_{\text{LL}}^{\text{sum}} = -\frac{s_0^g m_g^4}{g^2} f_1^{1-4\left(\frac{\gamma_0^g}{2b_0^g}\right)} = -\frac{s_0^g m_g^4}{g^2} f_1^{-1}, \\ f_1 = 1 + 2b_0^g g^2 \ln \frac{m_g}{M}, \quad (18) $$

where we also used Eq. (9). It is straightforward to check that Eq. (18) reproduces at all orders the coefficients in Eq. (12). The last equality in Eq. (18) tells that the LL series iterates simply like the (pure gauge) running coupling, but this is merely an accident of LO RG Eq. (9).

Similarly, after more algebra one can resum formally the NLL series. Adapting results from [31] we obtain

$$ P_{\text{NLL}}^{\text{sum}} = \frac{-s_0^g m_g^4}{g^2 f_2^{4A_0-1}} [R(f_2)]^B \left( 1 - \frac{a_{1,1} g^2}{s_0^g f_2} - \frac{a_{2,2} g^4}{s_0^g f_2^2} \right), \quad (19) \\ R(f_2) = [1 + g^2 b_1^g / (b_0^g f_2)] / (1 + g^2 b_1^g / b_0^g), \\ f_2 = 1 + [2b_0^g g^2 + 2(b_1^g - \gamma_0^g b_0^g) g^4] \ln \frac{m_g}{M} + \mathcal{O}(g^6), \quad (20) $$

and $A_0 = \gamma_0^g/(2b_0^g)$, $A_1 = \gamma_1^g/(2b_1^g)$, $B = 4(A_1 - A_0)$. The exact $f_2$ expression, reproducing Eq. (14) to all orders from Eq. (19), is given in [22]. Equation (20) gives numerically good approximations as long as the coupling is not too large ($\alpha_s \lesssim 0.5$).

A few remarks are worth making regarding the input content of Eq. (19): (i) it is rather generic, but numerically relies on the lowest order “purely soft” NLL coefficient $a_{2,1} \propto p_{-1}$ in Eq. (3), calculated in [19]. Accordingly, Eq. (19) does not include the (presently unknown) QCD mixed soft-hard NLL contributions; see Ref. [19]. (ii) One obtains $a'_{1,1} = a_{1,1} - s_1^g$ with $a_{1,1}$ defined in Eq. (5), since the NLO subtraction coefficient $s_1^g$ in Eq. (17) contributes a correction to $a_{1,1}$. (iii) $a_{2,2}$ in Eq. (19) incorporates the $\mathcal{O}(g^2 m_g^4)$ nonlogarithmic term $p_0$ in Eq. (3): the precise connection between the parameters in Eq. (19) and $p_{-1}$, $p_0$ in Eqs. (3) and (10) is $a_{2,1} = N_c/(4\pi)a_{1,0}(2p_{-1})$, $a_{2,2} = -N_c/(4\pi)a_{1,0}p_0$, after modifying $p_{-1}$, $p_0$ in Eq. (10).

**Soft and hard pressure matching.**—The massive RG construction above basically concerns the soft pressure contributions, thus $\ln m_g/M \to \ln m_E/M_s$, with overall factor $m_E^4 \sim \alpha_s^2$; see Eq. (2). While the hard contribution in Eq. (1), known exactly only at $\alpha_s^2$ order, is added perturbatively to Eqs. (18) and (19). The latter RGI expressions formally cancel the soft scale dependence up to neglected $\mathcal{O}(g^2 m_E^4)$, $\mathcal{O}(g^4 m_E^4)$ terms, respectively. At $\mathcal{O}(\alpha_s^2)$, it is equivalent to the complete soft scale cancellation operating in the factorization picture; thus, we can choose any $M_s$. While at $\alpha_s^{p \ge 3}$ orders, partial ignorance of hard and mixed contributions incites us to take $M_s \sim \mathcal{O}(m_E)$ in the $\ln m_E/M_s$ terms.

Another feature to account for while combining the soft and hard contributions is that Eq. (19) entails a nonlogarithmic term $\sim \alpha_s^2 a_{1,1}$ with $a_{1,1}$ from Eq. (17), obviously different from the genuine “soft +hard” $\alpha_s^2$ term in Eq. (1). Since only soft contributions are RG-resummed, to avoid wrong contaminations in nonlogarithmic hard terms from RG-induced NLL soft terms, one should perturbatively subtract $a_{1,1} m_E^4$ from the total expression [which does not affect higher order NLL terms generated by $a_{1,1}$ within Eq. (19)].

As a last important subtlety, note that although the subtraction terms as above determined, Eqs. (16) and (17), are sufficient to define LO and NLO RGI pressures
Eqs. (18) and (19), actually the complete integration of Eq. (15) entails an extra boundary condition (see Ref. [22]),

$$P_{\text{RGI}}^{\text{BC}} = m_g^4(M_0) \left( \frac{s_0^g}{g^2(M_0)} + s_1^g \right), \quad (21)$$

of similar form as the subtraction terms but involving a (boundary) scale $M_0 \neq M$. We can use this freedom to set $M_0$ such that Eq. (21) provides an appropriate EFT matching [27] of the soft pressure to the full one at $O(\alpha_s^2)$, i.e., with $M_0 \sim O(\mu)$, consistently also with the Stefan-Boltzmann limit.

Collecting all, our final pressure expression is obtained upon formally replacing $P_{\alpha_s^3}^{\text{soft}}$ in Eq. (1) by

$$P_{\alpha_s^{p \ge 3}}^{\text{sum}} = P_{(\text{N})\text{LL}}^{\text{sum}} + P_{\text{RGI}}^{\text{BC}} - m_E^4 \left( a_{1,1} + a_{1,0} \ln \frac{m_E}{M_s} \right) - P_{\alpha_s^3}^{\text{match}}, \quad (22)$$

with $a_{1,1} = 0$, $s_1^g \equiv 0$ for $P_{\text{LL}}^{\text{sum}}$, and $P_{\alpha_s^3}^{\text{match}}$ given in Eq. (B10) in [22]. The last three terms in Eq. (22) are required to match Eq. (1) consistently (i.e., without double counting), so that all RG-induced extra terms are $O(\alpha_s^3)$.

**Numerical results and comparisons.**—In Fig. 2, the RGI LL and NLL resummed pressures from Eq. (22) [32] are compared to the present state-of-the-art Eqs. (1), (3) as a function of $\mu_B = 3\mu$. The central scale values $M_h = 2\mu$ and the $\mu \le M_h \le 4\mu$ remnant scale dependence are illustrated for the different quantities, using in Eqs. (1) and (22) the exact NLO QCD running coupling $\alpha_s(M_h)$ with [33] $\Lambda_{\overline{\text{MS}}} \sim 0.32$ GeV. For sensible comparisons we also adopt the minimal sensitivity-determined [34] soft scale in [18], $M_s \sim 0.275 m_E$. Finally, we fix $M_0 = 2\mu$ in Eq. (21), a

<figure><img src="image_4.png" /><figcaption>FIG. 2. NNLO + soft N$^3$LO Eqs. (1) and (3) versus NNLO + soft N$^3$LO + RGI LL, NLL resummed pressures, as function of $\mu_B = 3\mu$, with $\mu \le M_h \le 4\mu$, $M_s \approx 0.275 m_E$, $M_0 = 2\mu$. For the RGI NLL resummed pressure, $M_s$ variations within [$M_s/2, 2M_s$] are shown in addition as darker bands.</figcaption></figure>

natural choice as it calibrates Eq. (22) to the central $M_h$ values of the NNLO pressure.

Note first importantly that the sole LL resummation, given by Eq. (18) with $-s_0^g \to a_{1,0}g^2 \ln m_g/M$, gives a sizeably reduced scale dependence compared to the NNLO pressure Eq. (1) at $O(\alpha_s^2)$, as Eq. (18) induces positive $\alpha_s^{p \ge 3}$ contributions partly canceling the negative $\alpha_s^2$ coefficient in Eq. (1). However, this effect is approximately cancelled once including $p_{-1}$, $p_0 \alpha_s^3$-order terms Eqs. (4) and (10). Next, for the LL and NLL RGI pressures, deviations from the state-of-the-art (“NNLO + soft N$^3$LO” in Fig. 2) are noticeable. The central scale ($M_h = 2\mu$) RGI pressure is slightly higher for fixed $\mu$ values, with very moderate differences between LL and NLL pressures. Importantly, the remnant scale dependencies of the resummed pressures are reduced as compared to NNLO + soft N$^3$LO results: only slightly for the LL pressure, due to cancellations with $p_{-1}$, $p_0 O(\alpha_s^3)$ terms, but significantly for the NLL one, both for $M_h$ and $M_s$ variations [35], due to $\alpha_s^{p \ge 3}$ terms induced both by NLL and the RGI-restoring terms in Eqs. (19) and (22).

In conclusion, we have obtained compact explicit expressions for the all order “double” resummations of LL and NLL soft contributions to the cold and dense QCD pressure that go well beyond previously established results. Our RG resummation construction moreover gives clearly improved residual scale dependence. This should provide improved control toward lower $\mu_B$ values to match with the extrapolated EOS from the nuclear matter density region. We thus anticipate that our present results may have strong implications once embedded within a more realistic EOS. It is not difficult to extend our framework to include different chemical potential and nonzero masses for the quarks in order to more realistically describe the EOS relevant for beta equilibrium and neutron star properties, which we leave for future investigation.

We thank Marcus B. Pinto and Aleksi Vuorinen for discussions.

[1] J. P. Blaizot, E. Iancu, and A. Rebhan, Thermodynamics of the high temperature quark gluon plasma, in *Quark Gluon Plasma* 3, edited by R. C. Hwa and X.-N. Wang (World Scientific, Singapore, 2004), pp. 60–122; U. Kraemmer and A. Rebhan, Rep. Prog. Phys. 67, 351 (2004).

[2] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Szabo, Nature (London) 443, 675 (2006); Y. Aoki, S. Borsanyi, S. Durr, Z. Fodor, S. D. Katz, S. Krieg, and K. K. Szabo, J. High Energy Phys. 06 (2009) 088; Szabolcs Borsányi, Zoltán Fodor, Christian Hoelbling, Sándor D. Katz, Stefan Krieg, Claudia Ratti, and Kálmán K. Szabó (Wuppertal-Budapest Collaboration), J. High Energy Phys. 09 (2010) 073; A. Bazavov, T. Bhattacharya, M. Cheng, C. DeTar, H. T. Ding, S. Gottlieb, R. Gupta, P. Hegde, U. M. Heller, F. Karsch et al., Phys. Rev. D 85, 054503 (2012).
[3] P. de Forcrand, Proc. Sci., LAT2009 (2009) 010 [arXiv:1005.0539]; G. Aarts, J. Phys. Conf. Ser. **706**, 022004 (2016).

[4] G. Baym, T. Hatsuda, T. Kojo, P. D. Powell, Y. Song, and T. Takatsuka, Rep. Prog. Phys. **81**, 056902 (2018).

[5] J. O. Andersen, E. Braaten, and M. Strickland, Phys. Rev. Lett. **83**, 2139 (1999); J. O. Andersen, E. Braaten, and M. Strickland, Phys. Rev. D **61**, 074016 (2000).

[6] J. O. Andersen, E. Braaten, E. Petitgirard, and M. Strickland, Phys. Rev. D **66**, 085016 (2002).

[7] J. O. Andersen, L. E. Leganger, M. Strickland, and N. Su, J. High Energy Phys. 08 (2011) 053; N. Haque, A. Bandyopadhyay, J. O. Andersen, M. G. Mustafa, M. Strickland, and N. Su, J. High Energy Phys. 05 (2014) 027.

[8] J. Ghiglieri, A. Kurkela, M. Strickland, and A. Vuorinen, Phys. Rep. **880**, 1 (2020).

[9] Y. Fujimoto and K. Fukushima, Phys. Rev. D **105**, 014025 (2022).

[10] J. L. Kneur, M. B. Pinto, and T. E. Restrepo, Phys. Rev. D **100**, 114006 (2019).

[11] J. L. Kneur, M. B. Pinto, and T. E. Restrepo, Phys. Rev. D **104**, 034003 (2021); **104**, L031502 (2021).

[12] B. A. Freedman and L. D. McLerran, Phys. Rev. D **16**, 1147 (1977); **16**, 1169 (1977).

[13] E. S. Fraga and P. Romatschke, Phys. Rev. D **71**, 105014 (2005).

[14] M. Laine and Y. Schröder, Phys. Rev. D **73**, 085009 (2006).

[15] A. Kurkela, P. Romatschke, and A. Vuorinen, Phys. Rev. D **81**, 105021 (2010).

[16] A. Vuorinen, Phys. Rev. D **68**, 054017 (2003).

[17] T. Gorda, A. Kurkela, P. Romatschke, M. Säppi, and A. Vuorinen, Phys. Rev. Lett. **121**, 202701 (2018).

[18] T. Gorda, A. Kurkela, R. Paatelainen, S. Säppi, and A. Vuorinen, Phys. Rev. Lett. **127**, 162003 (2021).

[19] T. Gorda, A. Kurkela, R. Paatelainen, S. Säppi, and A. Vuorinen, Phys. Rev. D **104**, 074015 (2021).

[20] E. Braaten and R. D. Pisarski, Phys. Rev. D **45**, R1827 (1992).

[21] S. Mogliacci, J. O. Andersen, M. Strickland, N. Su, and A. Vuorinen, J. High Energy Phys. 12 (2013) 055.

[22] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevLett.129.212001 for additional details.

[23] J. C. Collins, Renormalization (Cambridge University Press, Cambridge, England, 1984).

[24] J. L. Kneur and M. B. Pinto, Phys. Rev. Lett. **116**, 031601 (2016); Phys. Rev. D **92**, 116008 (2015).

[25] M. A. L. Capri, D. Dudal, J. A. Gracey, V. E. R. Lemes, R. F. Sobreiro, S. P. Sorella, and H. Verschelde, Phys. Rev. D **72**, 105016 (2005).

[26] M. A. L. Capri, D. Dudal, J. A. Gracey, V. E. R. Lemes, R. F. Sobreiro, S. P. Sorella, and H. Verschelde, Phys. Rev. D **74**, 045008 (2006).

[27] H. Georgi, Annu. Rev. Nucl. Part. Sci. **43**, 209 (1993); A. V. Manohar, Lectures at Les Houches summer school, arXiv:1804.05863.

[28] V. P. Spiridonov and K. G. Chetyrkin, Sov. J. Nucl. Phys. **47**, 522 (1988); see for a recent review P. A. Baikov and K. G. Chetyrkin, Proc. Sci., RADCOR2017 (2018) 025.

[29] B. M. Kastening, Phys. Rev. D **54**, 3965 (1996).

[30] Since $s_0^g, s_1^g$ are obtained from $a_{1,0}, a_{2,1}$, respectively, by applying the RG Eq. (6), we conveniently recast the LL, NLL series such that $s_0^g$ defines their first term, i.e., with $p \ge 1$ in Eq. (12) and $a_{0,0} \equiv -s_0^g$.

[31] J. L. Kneur, Phys. Rev. D **57**, 2785 (1998).

[32] See also Eqs. (B11),(C6) in [22] for compact expressions.

[33] M. Tanabashi et al. (Particle Data Group), Phys. Rev. D **98**, 030001 (2018).

[34] P. M. Stevenson, Phys. Rev. D **23**, 2916 (1981); Nucl. Phys. B**203**, 472 (1982).

[35] Note that $M_s \sim 0.275m_E$ maximizes the $M_h$ variations. As we also checked, uncertainties from varying $M_0$, for fixed $M_h$, are much smaller than the $M_h$ variations in Fig. 2.