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

HAL
open science

# Renormalization group improved pressure for cold and dense QCD

Jean-Loïc Kneur, Marcus Benghi Pinto, Tulio Eduardo Restrepo

► To cite this version:

Jean-Loïc Kneur, Marcus Benghi Pinto, Tulio Eduardo Restrepo. Renormalization group improved pressure for cold and dense QCD. Physical Review D, 2019, 100 (11), pp.114006. 10.1103/Phys-RevD.100.114006. hal-02302945

HAL Id: hal-02302945

https://hal.science/hal-02302945

Submitted on 1 Jun 2022

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Renormalization group improved pressure for cold and dense QCD

Jean-Loïc Kneur,¹,* Marcus Benghi Pinto,²,³,† and Tulio Eduardo Restrepo,²,‡

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

²Departamento de Física, Universidade Federal de Santa Catarina,
Florianópolis, Santa Catarina 88040-900, Brazil

³Department of Physics, University of Colorado, Boulder, Colorado 80309, USA

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

(Received 28 August 2019; published 2 December 2019)

We apply the renormalization group optimized perturbation theory (RGOPT) to evaluate the QCD (matter) pressure at the two-loop level considering three flavors of massless quarks in a dense and cold medium. Already at leading order ($\alpha_s^0$), which builds on the simple one loop (RG resummed) term, our technique provides a nontrivial nonperturbative approximation which is completely renormalization group invariant. At the next-to-leading order the comparison between the RGOPT and the perturbative QCD predictions shows that the former method provides results which are in better agreement with the state-of-the-art *higher order* perturbative results, which include a contribution of order $\alpha_s^3 \ln^2 \alpha_s$. At the same time one also observes that the RGOPT predictions are less sensitive to variations of the arbitrary $\overline{\text{MS}}$ renormalization scale than those obtained with perturbative QCD. These results indicate that the RGOPT provides an efficient resummation scheme which may be considered as an alternative to lattice simulations at high baryonic densities.

DOI: 10.1103/PhysRevD.100.114006

## I. INTRODUCTION

First principles evaluations aiming to describe the properties of strongly interacting matter at finite temperatures and/or baryonic densities are highly complicated by the inherent nonlinear and nonperturbative characteristics displayed by quantum chromodynamics (QCD). Nevertheless, at least in regimes of vanishing baryonic densities which concerns high energy heavy ion collisions, this fundamental theory can nowadays be successfully described by numerical lattice simulations (LQCD) [1]. However, the famous sign problem [2] for nonzero chemical potential is still preventing the method to be reliably applied to regimes of intermediate temperatures and baryonic densities which are relevant to experiments such as the *beam energy scan* (BES) at the Relativistic Heavy Ion Collider (RHIC) as well as CBM at FAIR and NICA at JINR which aim to locate the eventual QCD critical endpoint. At the same time the knowledge of an equation of state (EoS) that faithfully describes the cold and dense regime is necessary for an accurate description of

compact stellar objects. Unfortunately, for the reasons alluded above, LQCD cannot yet furnish such an EoS so that in general the problem is partially circumvented in different ways such as by using chiral effective theories (CET) [3] at low densities and perturbative QCD (pQCD) [4] at ultrahigh densities. As an alternative to these two (first principles) analytical approaches one may use effective quark models such as the MIT bag model [5], the Nambu–Jona-Lasinio model (NJL) [6], as well as the quark-meson model (QMM) [7] among others. In principle, pQCD applications should be carried out at extremely high densities where the asymptotic freedom property assures that the QCD coupling, $\alpha_s$, is small enough to justify the use of such an approximation. A seminal pQCD work by Freedman and McLerran [8,9] has provided the next-to-next-to-leading order (NNLO) pressure for massless quarks at vanishing temperatures and finite chemical potentials. The result has then been further refined so as to include thermal effects [10,11] and finite quark masses [12–14] apart from being rederived in a way compatible with the more modern $\overline{\text{MS}}$ renormalization scheme [15]. After more than four decades, a new perturbative order has recently been evaluated in Ref. [16], where the authors have determined the coefficient of the leading-logarithmic contribution at N3LO: $\alpha_s^3 \ln^2 \alpha_s$. Since the leading-logarithm soft contribution at N3LO evaluated in that work gives a negligible correction to the NNLO it was concluded that using pQCD result as an *ab-initio* input in calculations of the properties of

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³.

<footnote>jean-loic.kneur@umontpellier.fr</footnote>

<footnote>marcus.benghi@ufsc.br</footnote>

<footnote>tulio.restrepo@posgrad.ufsc.br</footnote>
neutron stars [13,17–20] as well as simulation gravitational-wave signals from neutron-star mergers is well justified. Nevertheless, for our present purposes it is also important to remark that in the evaluations performed in Ref. [16] the authors have chosen the $\overline{\text{MS}}$ arbitrary scale ($M$) to be $2\mu$ (where $\mu$ is the quark chemical potential). However, it is well known that physical observables evaluated with standard perturbation theory, as well as those obtained with resummation methods such as hard thermal loop perturbation theory (HTLpt) [21–23], can be very sensitive to renormalization scale variations. Moreover, it has been observed notably at finite temperature [4,23,24] that the latter scale sensitivity even *increases* when successive terms in the weak-coupling expansion are considered, which is an odd result as far as thermodynamical observables, such as the pressure, are concerned. At vanishing temperatures and finite baryonic densities, the scale dependence of the QCD pressure at NLO and NNLO has been explicitly investigated in Ref. [13]. The results show that the pressure has a rather large renormalization scale dependence, especially below the quark chemical potential $\mu \sim 1$ GeV, which corresponds to a baryon density $\sim 10^2 \rho_0$ where $\rho_0 \sim 0.16 \text{ fm}^{-3}$ represents the nuclear mass density. Such dependence indicates that the eventual nonperturbative effects remain quite important in the lower density range relevant to neutron stars. We remark that the renormalization scale dependence is even worse with the HTLpt resummation at finite $T$, where the calculations have been pushed to the three loop level, predicting results which agree with LQCD but only when the central scale $M = 2\pi T$ is used at $\mu = 0$ [22], while exhibiting a very large variation from $M = \pi T$ to $M = 4\pi T$. This is a clear indication that renormalization group (RG) properties have not been properly addressed in the perturbative and HTLpt resummed evaluations. The results displayed in Ref. [23] show that this unfortunate situation persists at finite densities and finite temperatures when the scale is varied around $M = 2\pi \sqrt{T^2 + \mu^2/\pi}$.

In the present work we consider an alternative resummation method which incorporates RG properties to evaluate the QCD pressure at $T = 0$ and finite $\mu$ values at the two loop level. This technique, which provides nonperturbative approximations, has been dubbed RGOPT (renormalization group optimized perturbation theory), and can be viewed as an extension of the standard optimized perturbation theory (OPT) [25,26] and the screened perturbation theory (SPT) [27] (both related to the so called linear $\delta$ expansion (LDE) [28]). Remark also that the HTLpt can be seen as the gauge-invariance compatible version of the OPT/SPT. Initially the RGOPT was employed [29] at vanishing temperatures and densities in the Gross-Neveu (GN) model [30]. Then, it has been applied to QCD at $T = \mu = 0$ to evaluate the basic scale $\Lambda_{\overline{\text{MS}}}$ [31,32] (equivalently the coupling $\alpha_s$), predicting values compatible with the world average [33]. The method has also been used to derive an accurate value of the quark

condensate [34]. More recently, it has been applied to the scalar $\lambda\phi^4$ theory [35,36], as well as to the nonlinear sigma model (NLSM) [37], producing results that show its compatibility with control parameters such as the temperature. The present paper is the first RGOPT application to (cold) in-medium QCD, for a nonzero chemical potential, so that one can analyze how the method performs in the regime of finite baryonic densities. The latter, despite being currently largely unaccessible to LQCD, is of utmost importance to the description of compact stellar objects. Our goal is twofold: first, we would like to check how our approach compares with the N3LO pQCD results recently obtained in Ref. [16]. Second, we aim to show how the scale dependence within the predicted pQCD pressure can be significantly reduced when the evaluations are performed within the RGOPT, a generic feature of the method. The paper is organized as follows. As a warm up, in the next section we review the RGOPT approach illustrating it with the $d=1+1$ massless Gross-Neveu model at $T = \mu = 0$. In Sec. III the method is used to evaluate the quark contribution to the QCD pressure at vanishing temperatures and finite densities up to the (RG optimized) NLO two-loop level. The optimization procedure and numerical results are presented in Sec. IV. Then in Sec. V we present our conclusions and perspectives.

## II. REVIEWING THE RGOPT WITH THE GN MODEL

The RGOPT belongs to a class of variational methods, reminiscent of the traditional Hartree approximation, which are particularly suitable to tackle infrared problems that plague massless theories. In this section the main steps of the approach will be recalled by performing a simple lowest order application to the massless Gross-Neveu model (GN) in two dimensions. More details and applications of the method can be found in Refs. [29,32,34–37]. The GN model is described by the Lagrangian density for a fermion field with $N$ components given by [30]

$$ \mathcal{L}_{\text{GN}} = \bar{\psi}(i\partial\!\!/\psi) + \frac{g_{\text{GN}}^2}{2} (\bar{\psi}\psi)^2. \quad (2.1) $$

The theory described by Eq. (2.1) is invariant under the transformation $\psi \to \gamma_5\psi$ displaying a discrete chiral symmetry (CS) in addition to having a global $O(2N)$ flavor symmetry. This simple renormalizable model has important common features with QCD, such as asymptotic freedom and a dynamically generated mass gap, among others. It is exactly solvable in the large-$N$ limit, and for arbitrary $N$ values the exact mass gap (at vanishing temperature) has been obtained [38] from Bethe ansatz methods. This allows to confront other nonperturbative approximation schemes that can include finite $N$ corrections (such as the RGOPT) to either the large- or finite-$N$ known results.
For convenience let us first rescale the four-fermion interaction as $g_{\text{GN}}^2 = g\pi/N$. To implement the RGOPT requires first to deform the interaction terms of Eq. (2.1) by introducing a Gaussian interpolating (mass) term and rescaling the coupling: in the case of a massless theory the RGOPT prescription is

$$ \mathcal{L}_{\text{GN}}^{\text{RGOPT}} = \mathcal{L}_{\text{GN}}(g \rightarrow \delta g) - m(1 - \delta)^a, \quad (2.2) $$

where $\delta$ is a book-keeping parameter interpolating between the free massive ($\delta = 0$) and interacting massless ($\delta = 1$) theory. Remark that setting $a = 1$ in Eq. (2.2) gives simply the “added and subtracted” variational mass prescription as adopted in the standard OPT/SPT/LDE [26–28]. In contrast a crucial feature of the RGOPT is to determine [29,31,32] the exponent $a$ from renormalization group consistency, giving generally $a \neq 1$, as we will recap below. Note that for the original massless model, the (free) propagator would normally be $S_F(p) = i(p)^{-1}$, while within the OPT or RGOPT approaches, any perturbative evaluations are first performed with a nonvanishing mass, thus providing an infrared regulator mass $m$, prior to the substitution Eq. (2.2) (the latter being most conveniently performed after a standard perturbative renormalization).

In the sequel of this section, to present a clearer overall picture of the approach we also restrict ourselves to the $T = \mu = 0$ case, since the main RGOPT features that we aim to recap are essentially determined by RG properties (thus by the renormalization aspects of the $T = \mu = 0$ part only). Once such RG properties are fixed, including the thermal and/or chemical potential contributions at a given order amounts to perform consistently the very same modifications as implied by Eq. (2.2) within those perturbatively calculated (massive) contributions.

We then start by evaluating the leading order $O(g^0)$ perturbative vacuum energy of the massive GN-model (more generally we could consider the pressure, with $T \neq 0$ and/or $\mu \neq 0$)

$$ \frac{\mathcal{E}_0^{PT}}{N} = i \int \frac{d^2 p}{(2\pi)^2} \ln(p^2 - m^2) + \mathcal{O}(g). \quad (2.3) $$

After renormalizing in the $\overline{\text{MS}}$-scheme (which at this lowest order amounts to simply a vacuum energy counterterm), one obtains

$$ \frac{\mathcal{E}^{PT}}{N} = -\frac{m^2}{2\pi} \left( \frac{1}{2} - L_m \right) + \mathcal{O}(g), \quad (2.4) $$

where $L_m = \ln(m/M)$ and $M$ is the arbitrary $\overline{\text{MS}}$ renormalization scale. Next consider the RG operator

$$ M \frac{d}{dM} = M \frac{\partial}{\partial M} + \beta(g) \frac{\partial}{\partial g} - \gamma_m(g) m \frac{\partial}{\partial m}, \quad (2.5) $$

with the normalization conventions for the RG coefficients [39]:

$$ \beta(g) = -b_0 g^2 - b_1 g^3 + \mathcal{O}(g^4), \quad (2.6) $$

$$ \gamma_m(g) = \gamma_0 g + \gamma_1^{\overline{\text{MS}}} g^2 + \mathcal{O}(g^3), \quad (2.7) $$

where $b_0 = 1 - 1/N$, $b_1 = -b_0/(2N)$, $\gamma_0 = 1 - 1/(2N)$, and $\gamma_1^{\overline{\text{MS}}} = -\gamma_0/(4N)$. The next step is to realize that Eq. (2.4) is not perturbatively RG-invariant: applying Eq. (2.5) to this expression gives a remnant term of leading order: $Md\mathcal{E}/dM = -m^2N/(2\pi) \neq O(g)$. However this rather well-known problem of a massive theory can be solved most conveniently by simply subtracting a (zero point) finite term in order to restore a RG invariant perturbative vacuum energy,<fnref n="1" /> that lead to the RG invariant (RGI) observable [32]

$$ \mathcal{E}^{\text{RGI}} = \mathcal{E}^{\text{PT}} - \frac{m^2}{g} s_0. \quad (2.8) $$

Now requiring Eq. (2.8) to satisfy the RG equations fixes the $s_0$ coefficient to

$$ s_0 = -N[2\pi(b_0 - 2\gamma_0)]^{-1} = \frac{N}{2\pi}. \quad (2.9) $$

The procedure is easily generalized most conveniently when taking higher perturbative order contributions into account by considering a perturbative subtraction

$$ -m^2 \sum_{k \ge 0} s_k g^{k-1}, \quad (2.10) $$

where the successive $s_i$ coefficients are fixed by requiring perturbative RG invariance, consistently including higher orders within the RG Eq. (2.5). This perturbative RG invariance restoration is of course not specific to the $d = 1 + 1$ GN model but more generic for any massive model, thus also in particular in four dimensions (see e.g., [40] for high order analysis in the $\phi^4$ theory). Now, incorporating those necessary subtraction terms, in order to start from a perturbatively RG invariant quantity, is an important necessary step prior to the specific RGOPT modification implied by Eq. (2.2). Next, performing the replacements Eq. (2.2) within a perturbative expression like (2.8) and doing a power expansion to order-$\delta^k$, one aims to recover formally the massless limit, $\delta \to 1$. But the latter (re) expansion leaves a remnant dependence on the (arbitrary) mass $m$ at any finite $\delta^k$ order, since the expression was initially perturbative. Indeed, applying the RGOPT replacements, at lowest $\delta^0$ order, to Eq. (2.8) gives

<footnote n="1">¹Alternatively one finds the very same results by requiring RG invariance at the level of bare expressions [29].</footnote>
$$ \frac{\mathcal{E}^{\text{RGOPT}}}{N} = -\frac{m^2}{2\pi} \left(\frac{1}{2} - L_m\right) - \frac{m^2}{Ng} (1 - 2a)s_0. \quad (2.11) $$

Now a crucial step is to realize that the resulting modified perturbative expression, Eq. (2.11), spoils the RG invariance in general, in particular for the simplest “added and subtracted mass” prescription $a = 1$, due to the drastic modification of the mass dependence. In contrast, the idea is to determine the interpolation exponent $a$ in Eq. (2.2) by requiring rather the reduced RG equation [29] to hold:

$$ f_{\text{RG}} \equiv M \frac{\partial \mathcal{E}^{\text{RGOPT}}}{\partial M} + \beta(g) \frac{\partial \mathcal{E}^{\text{RGOPT}}}{\partial g} \equiv 0, \quad (2.12) $$

in consistency with the massless limit being sought out. This uniquely fixes the exponent as

$$ a = \frac{\gamma_0}{b_0}, \quad (2.13) $$

a generic result also for other theories [32,34–37]. Moreover, the same value of $a$ is taken also when considering higher orders of the $\delta$-expansion, keeping the simple interpolating form of Eq. (2.2), since this exponent is universal (renormalization scheme independent).

Thus substituting $a = \gamma_0/b_0$ into Eq. (2.11) leads to the RGOPT lowest order result

$$ \frac{\mathcal{E}^{\text{RGOPT}}}{N} = -\frac{m^2}{2\pi} \left(\frac{1}{2} - L_m\right) + \frac{m^2}{2\pi g b_0}. \quad (2.14) $$

It is important to note that already at this lowest order the RGOPT-modified subtraction term clearly brings dynamical (RG) information through $g$ and $b_0$ apart from finite $N$ contributions (since $b_0 = 1 - 1/N$) to an otherwise trivial (free) vacuum energy. At this lowest order the final step consists in fixing the parameter $m$, still arbitrary at this stage, with an optimization prescription (MOP), similar to the so-called principle of minimal sensitivity (PMS) [25], defined by the stationarity condition

$$ f_{\text{MOP}} = \frac{\partial \mathcal{E}^{\text{RGOPT}}}{\partial m} \equiv 0 = \frac{m}{\pi} \left( \frac{1}{b_0 g} + L_m \right). \quad (2.15) $$

Apart from the trivial result $\bar{m} = 0$, one obtains

$$ \bar{m} = M \exp[-1/(b_0 g)], \quad (2.16) $$

which is clearly nonperturbative and explicitly RG invariant. Substituting $\bar{m}$ within Eq. (2.14), the vacuum energy is also RG invariant, and immediately gives the correct large-$N$ result, as was observed in Ref. [32]. To better appreciate these RGOPT features let us now compare this result with those obtained by the standard OPT/SPT as well as the large-$N$ approximations. As shown in Ref. [41], at order-$\delta^0$

the standard OPT/SPT vacuum energy (or equivalently pressure) has no information about the interactions since it is $g$-independent. Therefore the first non trivial contribution arises at next order-$\delta$ (two loop level) and by applying the MOP criterion one fixes the mass to

$$ \bar{m}_{\text{OPT}} = M \exp[-1/(g\gamma_0)], \quad (2.17) $$

which is not RG invariant. As for the $1/N$ expansion, the first nontrivial contribution appears at order-$N^0$ (the large-$N$ limit, LN) whose gap equation yields the well-known nonperturbative mass gap

$$ m^{\text{LN}} = M \exp(-1/g). \quad (2.18) $$

At this point a remarkable property of the RGOPT procedure over LN and standard OPT should be clear: it does produce a scale invariant nonperturbative result, which incorporates finite $N$ corrections, already at the one loop level. The same properties hold whenever adding in-medium contributions, because the latter are not affecting those RG properties which essentially rely on the vacuum contributions. Moreover, for $N \to \infty$ the RGOPT also reproduces the “exact” LN result, a consistency check of the reliability of the method. We point out that the latter property is also observed within the standard OPT since, as a particularity of the GN model, $\gamma_0 = b_0 \equiv 1$ at large-$N$.<fnref n="2" /> The LN limit is also reproduced when considering in-medium effects: for example for the $\phi^4$ model, quite remarkably the lowest order RGOPT pressure reproduces correctly [35] the (all-order) exact properties of the LN limit (that in the more standard large-$N$ derivation [43] involve the nontrivial resummation of “daisy” and “supers Daisy” graphs, associated with plasmon infrared divergent contributions). Yet for more involved theories such as QCD, one does not expect the lowest $\delta^0$-order RGOPT to be a very realistic approximation in general. This is because it essentially relies on lowest order RG quantities, while the other relevant contributions, e.g., in the pressure, are essentially those from a free theory at this order. Accordingly it appears sensible to go at least to the NLO order to get numerically more realistic results [32,34,36]. As we will illustrate in next sections this will be the case also for the QCD in-medium thermodynamic quantities considered in this work.

We will not proceed further with the GN model but to conclude this section we mention that the RGOPT recipe generalization to higher orders is rather straightforward, as will be better illustrated in the next sections with the in-medium QCD case. Once having implemented the relevant RG subtraction coefficients in Eq. (2.10), one performs the

<footnote n="2">²The fact that, for the GN as well as other theories, the OPT type of method does reproduce the $N \to \infty$ limit was observed long ago [42].</footnote>
RGOPT modification from Eq. (2.2) using the universal a value, Eq. (2.13), expanding this to $\delta^k$-order consistently with the original perturbative order considered, and taking the massless limit $\delta \to 1$. Finally one uses the RG Eq. (2.12) and (or) the MOP Eq. (2.15) to obtain “non-perturbative” approximations, in the sense that the resulting RG-consistent dressed mass is of order $\Lambda_{\overline{MS}}$ at $T = \mu = 0$ [32]. At nonvanishing temperatures the dressed mass also acquires a thermal dependence, but keeping its RG properties (see Refs. [35–37] for more detailed discussions). Ideally one would aim to solve the two Eqs. (2.15), (2.12) simultaneously to fix both a dressed running mass ($\bar{m}$) and a dressed running coupling ($\bar{g}$). However, as one proceeds to higher orders both equations often develop non linearities, so that an increasing number of solutions occur a priori, moreover not guaranteed to be all real-valued. These unwelcome features are indeed common with the other related OPT/SPT approaches. But in the RGOPT, Eq. (2.13) also crucially guarantees that the only acceptable solutions are those matching the standard perturbative behavior for $g \to 0$ at $T = 0$, a simple criteria that most often selects a unique solution, even at the highest (four-loop) order investigated so far [32,34]. Alternatively a less complete but often more handy RG compatible criterion requires to solving only the full RG Eq. (2.5), to fix the dressed mass $\bar{m}(g)$. Next the coupling (not yet fixed at this stage) is naturally traded for the ordinary running coupling at the relevant perturbative order, instead of being more nonperturbatively determined. Accordingly, the final physical quantities exhibit a more pronounced residual scale dependence, which can be interpreted as an estimate of the error introduced by this alternative procedure. Nevertheless, different applications have shown that this residual scale dependence is milder compared to the ones produced by standard PT and also by the related OPT/SPT approaches.

### III. RGOPT EVALUATION OF THE QCD QUARK PRESSURE

Let us now apply the RGOPT to the three flavor (dense matter) QCD up to order-$g$ (defining for convenience $g = 4\pi\alpha_s$), in the limit of vanishing temperatures and finite baryonic densities with $\mu_s = \mu_u = \mu_d \equiv \mu$, which is the

equilibrium condition for the massless case considered here. To thus treat properly the quark sector of QCD, the RGOPT requires to deforming the theory by rescaling the coupling (consistently for every standard QCD interaction terms) and a modified Gaussian interpolating (mass) term, following the prescription

$$ \mathcal{L}_{\text{QCD}}^{\text{RGOPT}} = \mathcal{L}_{\text{QCD}}|_{g \to \delta g} - m(1-\delta)^a \bar{\psi}_f \psi_f, \quad (3.1) $$

where $f = u, d, s$ is flavor index. The fermionic interpolating term proportional to $m$ is completely similar to the one previously discussed for the GN model. Note carefully that in order to compare with Ref. [16] in the present work we will investigate the case of vanishing current masses ($m_u = m_d = m_s = 0$), while $m$ in Eq. (3.2) above will become our variational mass upon implementing the RGOPT replacements, just as in the GN case illustrated in the previous Sec. II. (Accordingly $m$ represents a generic mass identical for the three flavors, in this initially SU(3) flavor symmetric approximation.) As a parenthetical remark, in principle a more complete and rather similar treatment of the gluon sector is possible, by following the hard thermal loop (HTL) prescription originally suggested by Braaten and Pisarski [44], that essentially introduces a gauge-invariant (nonlocal) effective Lagrangian properly describing a gluonic (thermal) “mass” term in the HTLpt approximation [21–23].

However, in the present work, which deals only with the $T = 0$ and $\mu \neq 0$ regime, we will apply the RGOPT to the quark sector only so that the gluon propagator, entering our evaluation at two-loop order, will be the usual (massless) one used in purely perturbative QCD (thus also with standard QCD interactions with quarks once the appropriate $\delta \to 1$ limit is taken, after the $\delta$-expansion following Eq. (3.1). This is justified by aiming to compare our results with the purely perturbative evaluation of the cold pressure such as in Ref. [16], also since the HTL-modified Lagrangian is supposed to play a crucial role more essentially once considering high temperature effects.

To order-$g$ the relevant contributions are displayed in Fig. 1. By combining the results of Ref. [34] for the vacuum ($\mu = 0$) contributions with those of Ref. [12] for the intermediate part one obtains the renormalized result

$$ P_{1,f}^{PT}(\mu) = -N_c \frac{m^4}{8\pi^2} \left(\frac{3}{4} - L_m\right) + \Theta(\mu^2 - m^2) \frac{N_c}{12\pi^2} \left[ \mu p_F \left(\mu^2 - \frac{5}{2}m^2\right) + \frac{3}{2}m^4 \ln\left(\frac{\mu + p_F}{m}\right) \right] \\ - \frac{d_A g}{4(2\pi)^4} m^4 \left(3L_m^2 - 4L_m + \frac{9}{4}\right) - \Theta(\mu^2 - m^2) \frac{d_A g}{4(2\pi)^4} \left\{ 3 \left[m^2 \ln\left(\frac{\mu + p_F}{m}\right) - \mu p_F\right]^2 - 2p_F^4 \right\} \\ - \Theta(\mu^2 - m^2) \frac{d_A g}{4(2\pi)^4} m^2 (4 - 6L_m) \left[ \mu p_F - m^2 \ln\left(\frac{\mu + p_F}{m}\right) \right], \quad (3.2) $$
<figure><img src="image_4.png" /><figcaption>FIG. 1. Feynman diagrams contributing to the perturbative quark pressure up to order-g.</figcaption></figure>

where $p_F = \sqrt{\mu^2 - m^2}$ is the Fermi momentum, $L_m = \ln(m/M)$, $d_A = N_c^2 - 1$, and $N_c = 3$. Now, to turn the above pressure into a RG invariant quantity, as explained in previous section, we subtract a finite “zero-point” contribution:

$$P_{1,f}^{PT}(\mu) \rightarrow P_{1,f}^{RGI}(\mu) = P_{1,f}^{PT}(\mu) - m^4 \sum_k s_k g^{k-1}. \quad (3.3)$$

Since our evaluations are being carried up to two-loop, order-g, it suffices to determine the first two coefficients $s_0$ and $s_1$ from applying the RG to the pressure (at $T = \mu = 0$), with the appropriate QCD $\beta$ and $\gamma_m$ RG functions. In our normalization conventions the QCD $\beta(g)$ and $\gamma_m(g)$ functions read

$$\beta(g \equiv 4\pi\alpha_s) = -2b_0g^2 - 2b_1g^3 + \mathcal{O}(g^4), \quad (3.4)$$

and

$$\gamma_m(g) = \gamma_0 g + \gamma_1 g^2 + \mathcal{O}(g^3), \quad (3.5)$$

where the coefficients are [45]

$$b_0 = \frac{1}{(4\pi)^2} \left( 11 - \frac{2}{3} N_f \right), \quad (3.6)$$

$$b_1 = \frac{1}{(4\pi)^4} \left( 102 - \frac{38}{3} N_f \right), \quad (3.7)$$

$$\gamma_0 = \frac{1}{2\pi^2} \quad (3.8)$$

and

$$\gamma_1^{\overline{\text{MS}}} = \frac{1}{8(2\pi)^4} \left( \frac{202}{3} - \frac{20}{9} N_f \right). \quad (3.9)$$

Applying the RG Eq. (2.5) to Eq. (3.2) and requiring the result to vanish up to higher $\mathcal{O}(g^2)$ terms determine the subtraction coefficients in Eq. (3.3) to

$$s_0 = -N_c [(4\pi)^2 (b_0 - 2\gamma_0)]^{-1}, \quad (3.10)$$

and

$$s_1 = -\frac{N_c}{4} \left[ \frac{b_1 - 2\gamma_1}{4(b_0 - 2\gamma_0)} - \frac{1}{12\pi^2} \right]. \quad (3.11)$$

Remark that the coefficients $s_k$, being determined solely from the vacuum contributions, do not depend on the mass nor on control parameters such as the temperature and chemical potential [35–37]. Next, to implement the actual RGOPT modification of interactions, we follow the substitution prescribed in Eq. (3.1). Like in the GN case the next step is to fix the exponent $a$, by expanding to leading order-$\delta^0$ and requiring the resulting pressure to satisfy the reduced RG Eq. (2.12), here applied to the QCD pressure. As expected this can be checked to yield the universal exponent:<fnref n="3" />

$$a = \frac{\gamma_0}{2b_0}, \quad (3.12)$$

in agreement with previous RGOPT applications to QCD [32,34]. The LO RGOPT pressure, per flavor, can then be written as

$$P_{0,f}^{\text{RGOPT}}(\mu) = N_c \frac{m^4}{(4\pi)^2 b_0 g} + \frac{N_c}{12\pi^2} \left[ \mu p_F \left( \mu^2 - \frac{5}{2} m^2 \right) + \frac{3}{2} m^4 \left( L_\mu - \frac{3}{4} \right) \right], \quad (3.13)$$

where $L_\mu \equiv \ln[(\mu + p_F)/M]$. Like for the GN model one can see that the RGOPT extra terms bring in information from RG dynamics through $g$ and $b_0$, to the otherwise trivial (free gas) perturbative result. Notice that in Eq. (3.13), assuming $\mu > m$ for most purpose below, the one-loop $\ln m$ terms of original Eq. (3.2) have canceled, as a result of considering the vacuum contributions given by the first term in Eq. (3.2), such that the $\ln(m/M)$ is consistently replaced by a $\ln[(\mu + p_F)/M]$ with the same coefficient.<fnref n="4" />

Next, considering the NLO RGOPT [i.e., taking $\delta \to 1$ in the $\delta^1$-order expansion of Eqs. (3.1), (3.3)] and after some algebra the modified pressure (per flavor) reads:

$$P_{1,f}^{\text{RGOPT}}(\mu) = P_{0,f}^{\text{RGOPT}}(\mu) - N_c \frac{m^4}{(4\pi)^2} \left(\frac{\gamma_0}{b_0}\right) \left(\frac{1}{b_0 g}\right) \\ + m^4 \left(2\frac{\gamma_0}{b_0} - 1\right) s_1 \\ + N_c \frac{m^2}{8\pi^2} \left(\frac{\gamma_0}{b_0}\right) [m^2(1 - 2L_\mu) + 2\mu p_F] \\ - \frac{gd_A}{4(2\pi)^4} \left[m^4\left(\frac{1}{4} - 4L_\mu + 3L_\mu^2\right) + \mu^2(\mu^2 + m^2) + m^2\mu p_F(4 - 6L_\mu)\right], \quad (3.14)$$

<footnote n="3">³Notice a trivial factor 1/2 difference as compared to Eq. (2.13) due to a convenient different normalization of $b_i$ in Eq. (3.4).</footnote>

<footnote n="4">⁴The very same cancellations occur in the original perturbative expression (3.2): this is not affected by RGOPT since the modification from (3.1) modifies all $\ln m$ terms similarly.</footnote>
where $P_{0,f}^{\text{RGOPT}}(\mu)$ is given by Eq. (3.13). Again, assuming $\mu > m$ for now on (except when explicitly mentioned below), we already simplified the $\ln m$ terms at two-loop order, as those originating from the vacuum contributions cancel exactly with those similar terms of the medium parts, so that Eq. (3.14) only depends on the combination<fnref n="5" /> $L_\mu = \ln[(\mu + p_F)/M]$. The LO and NLO RGOPT pressure are now ready to be optimized to generate nonperturbative approximation results as shown in the next section.

# IV. OPTIMIZATION PROCEDURE AND NUMERICAL RESULTS

## A. One-loop ($\delta^0$) RGOPT

Considering first the lowest ($\delta^0$) one-loop order result, let us recall that the constraint from the reduced RG Eq. (2.12) applied to the pressure, Eq. (3.13), has already been used to fix the exponent of the interpolating Lagrangian, see Eq. (3.12), such that by construction the pressure already satisfies $f_{\text{RG}} = 0$ exactly (at this order). Consequently the arbitrary mass $m$ may be fixed only by using the MOP optimization equation:

$$f_{\text{MOP}} = \frac{\partial P^{\text{RGOPT}}}{\partial m} \equiv 0. \qquad (4.1)$$

Considering first for simplicity solely the (one-loop) vacuum contribution in Eqs. (3.2) and (3.13) (reintroducing for this purpose consistently the $\ln m$ present at $\mu = 0$), one obtains

$$\bar{m}(\mu = 0) = \Lambda_{\overline{\text{MS}}} \sqrt{e}, \qquad (4.2)$$

where $\Lambda_{\overline{\text{MS}}} = Me^{-\frac{1}{2g\phi_0}}$ is the one loop QCD $\Lambda_{\overline{\text{MS}}}$ scale. Thus one obtains a nonperturbative mass, proportional to $\Lambda_{\overline{\text{MS}}}$, which is exactly RG invariant.

Including next the one-loop in-medium contribution from Eq. (3.13), we aim to use similarly the MOP Eq. (4.1) to now determine the $\mu$-dependent dressed mass $\bar{m}(\mu)$ (while the reduced RG equation is still automatically satisfied at this order for $\mu \neq 0$). At finite densities Eq. (4.1) is a little more involved due to the nonlinear $m$-dependence from $p_F(m) = \sqrt{\mu^2 - m^2}$ in Eq. (3.13). Yet, after simple algebra the formal solution may be cast into a compact form:

$$\bar{m}^2 = \mu^2 \left( \frac{\sqrt{1 + 4c(\bar{m}, \mu, g)} - 1}{2c(\bar{m}, \mu, g)} \right), \qquad (4.3)$$

where<fnref n="6" />

$$c(m, \mu, g) = \left( \frac{1}{2b_0 g} - \frac{1}{2} + L_\mu \right)^2. \qquad (4.4)$$

At this point we observe that the NLO subtraction $s_1$ in Eqs. (3.3), (3.11), while being strictly required for perturbative RG invariance only at two-loop order, is formally a one-loop $O(g^0)$ contribution. It appears thus sensible to include this known information from next RG order, which is straightforward and provides not surprisingly a somewhat more realistic one-loop improved approximation. Accordingly for $s_1 \neq 0$ the term $-1/2$ in Eq. (4.4) above is simply replaced by $-1/2 - 8\pi^2 s_1 = 11/84$ (for $n_f = 3$), which is the prescription used in the numerics below.

Equation (4.3) can be easily solved numerically but before doing that it is instructive to examine some of its properties in more detail. One can see first that the coupling $g \equiv g(M)$ and the renormalization scale $M$ only appear in the combination $1/(2b_0g) + L_\mu \approx 1/(2b_0g) + \ln(\mu/M) + \dots$, where the dots designate $M$-independent terms. Therefore, recalling that the (exact) one-loop running is defined as

$$g^{-1}(M) = g^{-1}(M_0) + 2b_0 \ln \left( \frac{M}{M_0} \right), \qquad (4.5)$$

for a reference scale $M_0$, it is immediate that Eq. (4.3) does not at all depend on $M$. Likewise it is easy to see that the (LO) RGOPT pressure Eq. (3.13) is itself exactly RG invariant at this one-loop order: formally replacing $m \to \bar{m}$ in its expression, $\bar{m}$ is RG invariant irrespectively of its numerical value, and the explicit $g(M)$ and $M$ in Eq. (3.13) only appear in the very same $M$-independent combination $1/(2b_0g(M)) + L_\mu$.

For small coupling, the optimum mass $\bar{m}$ admits a perturbative expansion $\bar{m}^2 \sim \mu^2(\text{constant} \times g + O(g^2))$, which has the expected form of an (in-medium) screening mass. Nevertheless, we insist at this point that $\bar{m}$ is not a physical mass, (and is not directly related to the Debye mass standard definition [46]), rather it represents an intermediate variational quantity whose sole purpose is to enter $P(\bar{m}, g, \mu)$, that defines the (optimized) physical pressure at a given order. In fact, except for very weak coupling, the first order expansion of $\bar{m}^2$ does not give a very good approximation of the exact $\bar{m}(\mu)$: indeed, instead of growing with no limits for arbitrary large coupling, as the purely perturbative approximation would naively suggest, the exact solution in Eq. (4.3) has the welcome property to be bounded, with $\bar{m}^2(g(M)) < \mu^2$ even for large $g(M)$.

<footnote n="5">⁵Those cancellations are however specific to the one- and two-loop level: at higher orders $\ln m$ and $\ln(\mu + p_F)$ appear independently, due to more “nested” divergences in the bare calculation [13].</footnote>

<footnote n="6">⁶Equation (4.3) suggests that $\bar{m}$ would be the solution of a simple quadratic equation if not for the nonlinear $m$ dependence entering $L_\mu = \ln[(\mu + p_F)/M]$. In Eq. (4.1) we have selected the solution $m^2 > 0$, while the other solution with $\sqrt{\cdots} \to -\sqrt{\cdots}$ is unphysical, giving always $m^2 < 0$.</footnote>
(therefore consistent with the basic assumptions of the in-medium contributions). The numerical solution $\bar{m}$ at LO RGOPT from Eq. (4.3) as a function of $\mu$ is illustrated in Fig. 2, which among other features evidently confirms its exact scale invariance properties.

## B. NLO two-loop ($\delta^1$) RGOPT: In-medium contribution

At NLO $O(g)$, it turns out that Eqs. (4.1) and (2.12) do not have real solutions for arbitrary chemical potential values. As already discussed above in Sec. II this is expected to happen in general, starting at NLO order, due to nonlinear dependences in the mass, if we insist to solve those equations exactly. Therefore, one could try less rigidly to solve the sole complete RG equation, Eq. (2.5), for $\bar{m}(g)$, taking then for $g$ more conservatively the purely perturbative two-loop running coupling. Unfortunately only nonreal solutions appear also in this case, if the equation is solved for exact $m(g)$. Nevertheless this situation can be remediated, at the price of introducing one extra parameter, to be fixed however by a well-defined prescription. Following Ref. [32], the idea is to modify the perturbative coefficients, expecting in this way to recover real solutions. But the modification should not be arbitrary, and should be RG compatible, so a presumably sensible prescription is to perform a (perturbative) renormalization scheme change (RSC). With a little insight, since one is mainly concerned with mass optimization, a simplest RSC can be defined by modifying only the mass parameter according to

$$m \rightarrow m'(1 + B_1 g + B_2 g^2 + \dots), \quad (4.6)$$

where the $B_i$ parametrize arbitrary scheme changes from the original $\overline{\text{MS}}$-scheme.<fnref n="7" /> For an exactly known function of $m$ and $g$ Eq. (4.6) would just be a change of variable not affecting physical results. While for a perturbative series truncated at order $g^k$, different schemes differ formally by remnant term of order $O(g^{k+1})$, such that the difference between two schemes is expected to decrease at higher orders for sufficiently weak coupling value. Now since we aim to solve optimization equations for “exact” $m$ and $g$ dependence, Eq. (4.6) actually modifies those purposefully, when now considered as constraints for the arbitrary mass $m'$. Furthermore we vary only one RSC parameter consistently at the same perturbative order, such that the relevant form of Eq. (4.6) is $m \rightarrow m'(1 + B_2 g^2)$: thus upon re-expanding to order-$g$ one can easily see that the net RSC modification to the pressure is to add the extra term $-4gm^4s_0B_2$ at two-loop order-$g$ (and simply renaming

$m' \rightarrow m$ afterwards the mass parameter to be determined to avoid excessive notation changes).

Clearly a definite prescription is required in order to fix $B_2$. Accordingly, one requires [32] the RSC to give the real $\bar{m}$ solution the closest to the original $\overline{\text{MS}}$-scheme: that is mathematically expressed by requiring the “contact” of the two curves (i.e., collinearity of the vectors tangent) parametrizing the relevant MOP and RG equations, considered as functions of $(m, g)$:

$$f_{\text{RSC}} = \frac{\partial f_{\text{RG}}}{\partial g} \frac{\partial f_{\text{MOP}}}{\partial m} - \frac{\partial f_{\text{RG}}}{\partial m} \frac{\partial f_{\text{MOP}}}{\partial g} \equiv 0, \quad (4.7)$$

where $f_{\text{MOP}}$ and $f_{\text{RG}}$ are given respectively by Eqs. (4.1) and (2.12) (applied here to the QCD pressure). Thus, Eq. (4.7) provides an extra constraint which completely fixes the additional RSC parameter $B_2$. Moreover, one expects the RSC to remain reasonably perturbative, i.e., $B_2$ to be moderate, which may be verified a posteriori by inspecting that $B_2g^2 \ll 1$.

As a first simple illustration, let us consider only the vacuum contribution at $\mu = 0$. Applying Eq. (4.7), in conjunction with the MOP Eq. (4.1) and taking for $g$ the two-loop perturbative Eq. (4.8), and for a typical value $M \approx 1$ GeV, one then obtains $\bar{B}_2 \approx -0.00224$, $\ln(\bar{m}/M) \approx -0.331$, giving $\bar{m} \approx 700$ MeV $\approx 2.1\Lambda_{\overline{\text{MS}}}$, which may be compared to the LO RGOPT result Eq. (4.2). Note that $\alpha_s(M = 1$ GeV) $\approx 0.42$, but $|\bar{B}_2g^2| \sim 0.06$ is a very moderate deviation from original $\overline{\text{MS}}$-scheme.

We can now numerically optimize the NLO RGOPT pressure Eq. (3.14) including the in-medium contribution with $\mu \neq 0$, adopting the RSC additional prescription to recover real $\bar{m}(\mu)$ solutions for arbitrary $\mu$ values. Actually, rather than solving the full RG, as a numerically simpler variant we solve the MOP equation (4.3) at two-loop order for $\bar{m}(g)$, taking then for $g$ the purely perturbative running coupling, together with the RSC equation (4.7) to fix $\bar{B}_2(\mu)$ at NLO. (Clearly the optimal RSC parameter $\bar{B}_2$ will now be a nontrivial function of the chemical potential $\mu$, consistently determined by the optimization procedure). At two-loop the analytical expression of the MOP, Eq. (4.1), is more involved than its one-loop analogue, so that the algebra leading to the solution Eq. (4.3) is not as easy but it can be readily solved numerically.

In order to compare with the results given in Refs. [13,16] we will consider the scale variation $\mu \le M \le 4\mu$ besides the “central” scale $M = 2\mu$. The exact two-loop (2L) running coupling, analogue of the one-loop Eq. (4.5), is obtained by solving for $g(M)$ the relation

$$\ln \frac{M}{\Lambda_{\overline{\text{MS}}}} = \frac{1}{2b_0 g} + \frac{b_1}{2b_0^2} \ln \left( \frac{b_0 g}{1 + \frac{b_1}{b_0} g} \right), \quad (4.8)$$

for a given $\Lambda_{\overline{\text{MS}}}$ value (this also defines the (two-loop order) $\Lambda_{\overline{\text{MS}}}$ in our normalization conventions). Equivalently to

<footnote n="7">⁷Equation (4.6) has also the welcome property that it does not affect the definition of the reference QCD scale $\Lambda_{\overline{\text{MS}}}$, in contrast with a similar perturbative change on the coupling, see [32] for details.</footnote>
giving a $\Lambda_{\overline{\text{MS}}}$ value one can give a $g(M_0)$ at some reference scale, $M_0$. Here, we have chosen $\alpha_s(M_0 = 1.5 \text{ GeV}) = 0.326$ to compare precisely with the values adopted in Ref. [16]: with (4.8) this corresponds to $\Lambda_{\overline{\text{MS}}} \approx 0.335 \text{ GeV}$, a value indeed very close to the present world average [33].

The results for the optimized NLO mass $\bar{m}$ as functions of $\mu$ and for different renormalization scale choices are shown in Fig. 2 where they are also compared with the LO one-loop $\bar{m}$ value. One can see that, as already explained above, $\bar{m}(\mu)$ is exactly RG invariant at LO, because it only involves the scale invariant combination $[2b_0g(M)]^{-1} + \ln(\mu/M)$. In contrast the NLO $\bar{m}(\mu)$ definitely displays a residual scale dependence: even for the exact two-loop running, Eq. (4.8), the latter is not very surprisingly no longer exactly “matched” by the optimized NLO RGOPT mass. We will illustrate below that the RGOPT pressure, which represents the actual physical observable, shows a more moderate residual scale dependence. More generally the RGOPT construction only guarantees that the optimization does not spoil the perturbative RG invariance of the physical quantity considered, that means up to remnant scale-dependent terms of higher order $O(g^{k+1})$, if the original perturbative expression is available at order $g^k$.

Figure 3 illustrates the corresponding values of the RSC parameter combination $B_2(\mu)g^2(\mu)$, thus quantifying the departure from $\overline{\text{MS}}$-scheme. One can see that RSC remains reasonably perturbative, although the value of $|B_2(\mu)|$ needed to recover real solutions are increasing rapidly for smaller $\mu$ values for the lower renormalization scale $M = \mu$ (not surprisingly since in this region the running coupling $g(M)$ becomes dangerously large).

One is now in position to compute thermodynamical observables, such as the pressure and the quark number density which here will be respectively normalized by the

<figure><img src="image_5.png" /><figcaption>FIG. 2. The optimized mass as a function of the chemical potential at $M = \mu$, $2\mu$ and $4\mu$ at order-$g^0$ (dot-dashed) and order-$g$ (continuous). For the latter the upper curve corresponds to $M = \mu$, the central curve to $M = 2\mu$, and the lower curve to $M = 4\mu$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 3. The optimized RSC $\bar{B}_2(\mu)g^2(\mu)$ quantity as a function of the chemical potential at $M = \mu$, $2\mu$ and $4\mu$ at order-$g$.</figcaption></figure>

equivalent massless free gas quantities $P_{fg}$ and $\rho_{fg}$. These quantities, per flavor, are respectively

$$P_{fg} = N_c \frac{\mu^4}{12\pi^2}, \qquad (4.9)$$

and

$$\rho_{fg} = N_c \frac{\mu^3}{3\pi^2}. \qquad (4.10)$$

Let us then compare the $O(g^0)$ and $O(g)$ RGOPT results with the pQCD predictions at $O(g)$, $O(g^2)$ [13], as well as the most recent $O(g^3 \ln^2 g)$ [16]. For completeness, we recall that the relevant pQCD expression is [16]

$$\frac{P_{\text{pQCD}}}{P_{fg}} = 1 - \frac{2}{\pi}\alpha_S(M) - \alpha_S^2(M) \left\{ \begin{aligned} & 0.303964 \ln \alpha_S(M) \\ & + \left[ 0.874355 + 0.911891 \ln \left( \frac{M}{\mu} \right) \right] \\ & - 0.266075 \alpha_S^3(M) \ln^2 \alpha_S. \end{aligned} \right\} \qquad (4.11)$$

In Fig. 4 we show the normalized pressure predicted by the different order approximations at the central scale choice $M = 2\mu$ as adopted in Ref. [16]. The first thing to remark is that the RGOPT produces a nontrivial result already at order-$\delta^0 g^0$, but converging quite slowly to the free gas result as the quark chemical potential increases. Nevertheless this can already be seen as an improvement since at this same order $g^0$ the pQCD result for the normalized pressure would trivially be equal to the unity, i.e., the free gas limit. In fact the lowest order RGOPT cannot be expected to be a very realistic approximation in general, because it only relies on lowest order RG quantities, while the pressure dependence is essentially like the free gas one. The efficient resummation properties of RGOPT become more evident when one compares its result at NLO, order-$g$, with the pQCD ones at the same NLO order, since the figure shows that the NLO RGOPT
<figure><img src="image_7.png" /><figcaption>FIG. 4. The normalized pressure as a function of the chemical potential at the central scale $M = 2\mu$. pQCD results at orders $g$, $g^2$, and $g^3$ (LL term) are compared with the RGOPT results at orders $g^0$ (one loop) and $g$ (two loop).</figcaption></figure>

pressure actually appears in much better agreement with the higher order perturbative $g^2$ and $g^3 \ln^2 g$ predictions. Next we also analyze how the different approximations perform when the arbitrary renormalization scale is varied in the range $\mu \le M \le 4\mu$, as in Ref. [13] where the scale dependence of pQCD results at orders $g$ and $g^2$ have been analyzed.<fnref n="8" /> The results are compared in Fig. 5, where the RGOPT appears to moderately improve the scale uncertainty, at least in the range $\mu \gtrsim 1$ GeV, as compared with the same perturbative order $g$.

To assess more precisely the remnant scale dependence we plot in Fig. 6 the difference of the (normalized) pressures $\Delta P/P_{fg} = (P(M=4\mu) - P(M=\mu))/P_{fg}$ as function of $\mu$, for the three approximations illustrated in Fig. 5. The NLO RGOPT remnant scale dependence is moderately but clearly improved as compared to NLO pQCD for $\mu \gtrsim 0.9$ GeV (giving ~25% improvement e.g., for $\mu \approx 2$ GeV), while the NLO pQCD scale dependence appears somewhat smaller in the lower $\mu$ range $0.5 \lesssim \mu \lesssim 0.9$ GeV. Notice also that the NNLO pQCD pressure has a smaller scale dependence than NLO pQCD in a narrower and more perturbative range $\mu \gtrsim 1.5$ GeV. In contrast the RGOPT scale uncertainty is clearly better than the NNLO pQCD one in the full relevant $\mu$ range. We remark however that the smaller remnant dependence of NLO pQCD within the low-$\mu$ window ($0.5 \lesssim \mu \lesssim 0.9$ GeV) is merely a side effect of the NLO pQCD pressure dropping toward zero at lower $\mu$ values than the two other approximations, as is clear from Fig. 5. Indeed not surprisingly all three approximations exhibit a rapidly growing scale

<figure><img src="image_8.png" /><figcaption>FIG. 5. The normalized pressure as a function of the chemical potential. pQCD results at NLO order-$g$ and NNLO order-$g^2$ are compared with RGOPT at NLO order-$g$. In each case the upper curve corresponds to $M = 4\mu$ and the lower curve to $M = \mu$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 6. The remnant scale dependences defined by the differences $\Delta P/P_{fg} = (P(M=4\mu) - P(M=\mu))/P_{fg}$ of (normalized) pressures, as functions of the chemical potential $\mu$. pQCD results at NLO order-$g$ and NNLO order-$g^2$ are compared with RGOPT at NLO order-$g$.</figcaption></figure>

dependence for $\mu$ values approaching the region where $P(M \approx \mu)$ rapidly drops toward zero.<fnref n="9" /> But Fig. 6 also shows that the maximal remnant dependence reached at the respective $\mu_{min}$ values is smaller for the RGOPT than for NLO and NNLO pQCD. In any case one should keep in mind that, due to the adopted common renormalization scale choice $g(M=O(\mu))$, none of the approximations is much reliable in the nonperturbative region where $P/P_{fg}(M \approx \mu) \ll 1$ due to large coupling [note, e.g., that $\mu < 0.8$ already corresponds to $\alpha_s(M=\mu) > 0.5$ using Eq. (4.8)]. We thus conclude that, within the $\mu$ range where all the approximations are very reliable perturbatively, the NLO RGOPT remnant scale uncertainty is moderately but

<footnote n="8">⁸In the original study [13] a quite common approximate form of (4.8) was rather used, truncating terms beyond $O(\ln L/L^2)$, with $L \equiv \ln(M^2/\Lambda_{\text{MS}}^2)$. Here we compare the scale dependence by adopting the same exact two-loop running coupling (4.8) for all approximations, that tends to very slightly decrease the remnant scale uncertainty for all cases.</footnote>

<footnote n="9">⁹In Fig. 6 the three curves consistently start at their respective minimal $\mu_{min}$ values, defined such that $P(M \ge \mu_{min}) \ge 0$, compare with Fig. 5.</footnote>
clearly improved in relative comparison to both NLO and NNLO pQCD (considering also that standard pQCD at $T = 0$, $\mu \neq 0$ has anyway less severe scale dependence issues than in the high $T$ regime).

In principle, we could also include in our NLO RGOPT analysis the NNLO $gm^4s_2$ subtraction term of Eq. (3.3), since being formally of order-$g$, similarly to what was done at LO RGOPT [see the discussion after Eq. (4.4)]. The $s_2$ expression is available from [34] and clearly incorporates additional RG dependence from next (three-loop) RG order. However we have checked that considering $s_2 \neq 0$ at NLO scarcely changes our results (in contrast with the LO pressure where $s_1 \neq 0$ has a sizeable impact). In particular the scale dependence is not visibly affected, which signals that a reasonable stability has been reached at NLO order.

In Ref. [13] the authors also analyzed the predictions for the quark number density:

$$ \rho(\mu) \equiv \frac{dP(\mu)}{d\mu}, \qquad (4.12) $$

up to NNLO pQCD. Their results are reproduced and compared with our RGOPT predictions in Fig. 7. As in the case of the pressure a noticeable (but moderate) decrease of the scale “uncertainty” band occurs for the NLO RGOPT density (while the LO RGOPT results are again exactly RG invariant for the same reasons than the LO pressure). However the RGOPT scale dependence improvement is less pronounced than for the pressure, which can be traced to our use of the standard running coupling (having renounced, as explained above, to the more complete optimization of $g$ and $m$, due to non real and involved solutions). Indeed, for dense matter the imperfectly balanced scale dependence, from the contribution of the running $g(M)$, tends to be enhanced as compared to the pressure since taking $g(M \sim \mu)$ to obtain $\rho(\mu)$ in Eq. (4.12)

<figure><img src="image_10.png" /><figcaption>FIG. 7. The quark number density as a function of the chemical potential. pQCD results at NLO order-$g$ and NNLO order-$g^2$ are compared with the RGOPT at NLO order-$g$. In each case the upper curve corresponds to $M = 4\mu$ and the lower curve to $M = \mu$.</figcaption></figure>

involves a contribution $\propto \partial_M g(M)$ on top of the explicit derivative $\partial_\mu P$ term.

### C. A simpler alternative NLO RGOPT prescription

While the results in Figs. 4 and 5 clearly show a better agreement of NLO RGOPT with the state-of-the-art perturbative results, it may be regarded rather unsatisfactory to have to deal with the somewhat more involved RGOPT NLO prescription, that implies the additional constraint from RSC Eq. (4.7) to be altogether numerically solved to restore real solutions. Could we find a simpler and more transparent prescription, while still capturing the main features of the RGOPT approach? Indeed, a much simpler alternative that surely recovers a real $\bar{m}$ solution is simply to renounce to solving the RG or MOP equations exactly, by approximating the latter in a more perturbative fashion. (This, however, certainly loses a part of the resummation properties embedded in the “exact” solution, such that a slight degradation of the remnant scale dependence is to be expected). To explore this alternative we consider the full RG equation (2.5), in order to incorporate the most complete and consistent NLO RG content, but we approximate crudely its solution to its first perturbative (re) expansion order. Similarly as in the LO case, noting first that Eq. (2.5) would give a simply quadratic equation for $m^2$ in absence of the extra nonlinear $m$-dependence from $p_F$, this perturbative solution is simple:<fnref n="10" />

$$ \begin{aligned} \bar{m}^2 &= \frac{9}{7\pi^2} \mu p_F \left( 1 - \sqrt{1 - \frac{21 \mu^2}{64 p_F^2}} \right) g + \mathcal{O}(g^2) \\ &= \frac{9}{7\pi^2} \left( 1 - \frac{\sqrt{43}}{8} \right) g\mu^2 + \mathcal{O}(g^2). \end{aligned} \qquad (4.13) $$

Now, inserting this $\bar{m}$ expression into the NLO RGOPT pressure expression Eq. (3.14), with the running $g \to g(M)$ as previously from Eq. (4.8), gives the results shown in Fig. 8, which are compared with pQCD at NNLO including the four-loop (LL) results from Eq. (4.11) (originally obtained in Ref. [16]). We illustrate also the scale dependence for the range $\mu \le M \le 4\mu$ for the two expressions. One sees the quite remarkable agreement for the central scale choice $M = 2\mu$ (more precisely with less than ~1.5% differences for any $\mu > 0.6$ GeV), while the RGOPT scale “uncertainty” range is still slightly better even for this rather crude approximation. Concerning the scale dependence

<footnote n="10">¹⁰In Eq. (4.13) the factors 9/7, 21/64, $\sqrt{43}$ are simply $n_f = 3$ values of the specific combinations of RG coefficients $b_i$, $\gamma_i$ appearing in this expression. The explicitly scale-dependent term $\ln \mu/M$ only appears at next $g^2$-order. Note also that we eliminated the other solution, with $+\sqrt{\cdots}$, as it violates the necessary consistency $\bar{m}^2 \le \mu^2$ even for moderate $g$, and also does not fulfill the perturbative matching, $\ln \mu/\bar{m} \sim 1/(2b_0g)$, for $\mu \gg \bar{m}$.</footnote>
<figure><img src="image_11.png" /><figcaption>FIG. 8. The normalized pressure as a function of the chemical potential. pQCD results from Eq. (4.11) at NNNLO including the contribution $g^3 \ln^2 g$ (blue dotted curves) compared with the simpler alternative NLO RGOPT (black continuous curve). In each case the upper curve corresponds to $M = 4\mu$, the central curve to $M = 2\mu$, and the lower curve to $M = \mu$.</figcaption></figure>

band of pQCD including the highest available perturbative order result, Eq. (4.11), it hardly displays a visible difference with the sole NNLO, order-$g^2$, perturbative pressure as studied in Ref. [13]. But the net effect of the highest order last term of Eq. (4.11), being negative, is to shift down (very slightly) the values of the pressure for given $\mu$ and $M$ values. Examining Fig. 4 we further observe that going from NLO to NNLO pQCD there is a more pronounced decrease of the pressure for given $\mu$ values [which is clear from the globally negative NNLO $O(g^2)$ terms in Eq. (4.11)]. Now, in Fig. 4 the *exact* NLO RGOPT pressure values are sensibly lower than the other approximations, while in contrast the approximate NLO RGOPT pressure, obtained with the perturbative $\bar{m}$ Eq. (4.13), agrees quite neatly with Eq. (4.11). Accordingly one may hint from those comparisons that the “exact” NLO RGOPT result may be a more precise approximation than Eq. (4.11) to the even higher order perturbative pressure values.

V. CONCLUSIONS

In this work we have performed the first application of the RGOPT resummation to QCD when a control parameter, such as the chemical potential, is present. As discussed this technique generates nonperturbative approximations with consistent RG properties in a region of the QCD phase diagram which is currently unavailable to LQCD simulations. Our results have been compared to the state-of-the-art pQCD predictions that include a $\alpha_s^3 \ln^2 \alpha_s$ contribution. We have confirmed in this in-medium application the generic property that at lowest one-loop order this technique already captures nontrivial and RG invariant results for the pressure and the quark number density. Although numerically these lowest order results are a poor approximation in general, and converge quite slowly to the free gas result as $\mu$ increases, they exhibit the more efficient RGOPT

resummation since, at this same order, the pQCD prediction is trivial. At NLO order-$g$ (two loop level) and $M = 2\mu$ the RGOPT results appear to be a very good approximation as they show a much better agreement with the perturbative higher orders $O(\alpha_s^3 \ln^2 \alpha_s)$ than pQCD at the same order. Scale variations in the range $M = \mu - 4\mu$ also show that the method reduces the scale uncertainties (although moderately at two-loop order) as compared to pQCD, which is important as far as EoS suitable to describe neutron stars are concerned.

The scale uncertainty improvement from RGOPT thus appears less spectacular than for other models explored at two-loop orders at $T \neq 0, \mu = 0$ compared with standard perturbation and HTLpt [35–37]. But this is merely due to the fact that standard pQCD at $T = 0, \mu \neq 0$ has less severe remnant scale dependence issues (as already noted in Ref. [13]) than most other models have in the high $T$ regime. In contrast the NLO RGOPT scale uncertainty appears more similarly moderate in both regimes. As discussed in the text and in other applications (see, e.g., Ref. [35]) the appearance of a residual (mild in most cases) scale dependence is unavoidable within the RGOPT beyond LO. But it is also clear [35] that since RGOPT maintains by construction the most possible of (perturbative) RG invariance, generically the scale uncertainty bands observed at NLO should further shrink by considering the NNLO, $O(g^2)$, which should also provide a priori more accurate numerical results. We remark that by combining the three loop vacuum contributions of Ref. [34] with the in-medium contributions of Ref. [13] this is a feasible, although technically more involved analysis (regarding optimization), that we intend to address in a future investigation. Regarding the present application, where only massless quarks have been considered, our results indicate that this RG-consistent resummation method is suitable to treat dense and cold QCD. Note also that it can be easily extended to determine more realistic EoS (e.g., including massive quarks) which aim to describe neutron stars. Finally, the RGOPT interpolation should be extended to the gluonic sector for a more complete description especially when considering high temperature effects [47].

ACKNOWLEDGMENTS

This work was financed in part by Coordenação de Aperfeiçoamento de Pessoal de Nível Superior—(CAPES-Brazil)—Finance Code 001 and by INCT-FNA (Process No. 464898/2014-5). T. E. R. acknowledges Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq-Brazil) and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES-Brazil) for Ph.D. grants at different periods of time. M. B. P. and J.-L. K. thank Paul Romatschke for discussions and the Department of Physics at UC Boulder, where this work was completed, for the hospitality. Instituto Nacional de Ciência e Tecnologia de Física Nuclear e Aplicações
[1] Y. Aoki, G. Endrődi, Z. Fodor, S. D. Katz, and K. K. Szabó, Nature (London) **443**, 675 (2006); Y. Aoki, S. Borsányi, S. Dür, Z. Fodor, S. D. Katz, S. Krieg, and K. K. Szabó, J. High Energy Phys. 06 (2009) 088; S. Borsányi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabó (Wuppertal-Budapest Collaboration), J. High Energy Phys. 09 (2010) 073; A. Bazavov et al. (HotQCD Collaboration), Phys. Rev. D **85**, 054503 (2012).

[2] P. de Forcrand, Proc. Sci., LAT2009 (2009) 010; G. Aarts, J. Phys. Conf. Ser. **706**, 022004 (2016).

[3] R. Machleidt and D. Entem, Phys. Rep. **503**, 1 (2011).

[4] U. Kraemmer and A. Rebhan, Rep. Prog. Phys. **67**, 351 (2004).

[5] A. Chodos, R. L. Jaffe, K. Johnson, and C. B. Thorn, Phys. Rev. D **10**, 2599 (1974).

[6] Y. Nambu and G. Jona-Lasinio, Phys. Rev. **122**, 345 (1961); **124**, 246 (1961); M. Buballa, Phys. Rep. **407**, 205 (2005).

[7] M. Gell-Mann and M. Lévy, Nuovo Cimento (1955-1965) **16**, 705 (1960); R.-A. Tripolt, B.-J. Schaefer, L. von Smekal, and J. Wambach, Phys. Rev. D **97**, 034022 (2018).

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

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

[10] A. Ipp, K. Kajantie, A. Rebhan, and A. Vuorinen, Phys. Rev. D **74**, 045016 (2006).

[11] A. Kurkela and A. Vuorinen, Phys. Rev. Lett. **117**, 042501 (2016).

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

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

[14] E. S. Fraga, A. Kurkela, and A. Vuorinen, Astrophys. J. **781**, L25 (2014).

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

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

[17] A. Kurkela, E. S. Fraga, J. Schaffner-Bielich, and A. Vuorinen, Astrophys. J. **789**, 127 (2014).

[18] T. Gorda, Astrophys. J. **832**, 28 (2016).

[19] E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, Phys. Rev. Lett. **120**, 172703 (2018).

[20] E. R. Most, L. R. Weih, L. Rezzolla, and J. Schaffner-Bielich, Phys. Rev. Lett. **120**, 261103 (2018).

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

[22] J. O. Andersen, M. Strickland, and N. Su, Phys. Rev. Lett. **104**, 122003 (2010); J. High Energy Phys. 08 (2010) 113; J. O. Andersen, L. E. Leganger, M. Strickland, and N. Su, J. High Energy Phys. 08 (2011) 053.

[23] S. Mogliacci, J. O. Andersen, M. Strickland, N. Su, and A. Vuorinen, J. High Energy Phys. 12 (2013) 055; N. Haque, J. O. Andersen, M. G. Mustafa, M. Strickland, and N. Su, Phys. Rev. D **89**, 061701 (2014); N. Haque, A. Bandyopadhyay, J. O. Andersen, M. G. Mustafa, M. Strickland, and N. Su, J. High Energy Phys. 05 (2014) 027.

[24] P. Blaizot, E. Iancu, and A. Rebhan, in *Thermodynamics of the High Temperature Quark Gluon Plasma*, edited by R. C. Hwa et al., Advanced Series on Directions in High Energy Physics Vol. 6 (World Scientific, Singapore, 2003), p. 60.

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

[26] S. Chiku and T. Hatsuda, Phys. Rev. D **58**, 076001 (1998).

[27] F. Karsch, A. Patks, and P. Petreczky, Phys. Lett. B **401**, 69 (1997); J. O. Andersen, E. Braaten, and M. Strickland, Phys. Rev. D **63**, 105008 (2001); J. O. Andersen and M. Strickland, Phys. Rev. D **64**, 105012 (2001).

[28] V. I. Yukalov, Theor. Math. Phys. **28**, 652 (1976); W. E. Caswell, Ann. Phys. (N.Y.) **123**, 153 (1979); I. Halliday and P. Suranyi, Phys. Lett. B **85**, 421 (1979); R. P. Feynman and H. Kleinert, Phys. Rev. A **34**, 5080 (1986); H. Jones and M. Moshe, Phys. Lett. B **234**, 492 (1990); A. Neveu, Nucl. Phys. B, Proc. Suppl. **18**, 242 (1991); V. I. Yukalov, J. Math. Phys. (N.Y.) **32**, 1235 (1991); C. M. Bender, F. Cooper, K. A. Milton, M. Moshe, S. S. Pinsky, and L. M. Simmons, Phys. Rev. D **45**, 1248 (1992); H. Yamada, Z. Phys. C **59**, 67 (1993); A. Sissakian, I. Solovtsov, and O. Solovtsova, Phys. Lett. B **321**, 381 (1994); H. Kleinert, Phys. Rev. D **57**, 2264 (1998); Phys. Lett. B **434**, 74 (1998).

[29] J.-L. Kneur and A. Neveu, Phys. Rev. D **81**, 125012 (2010).

[30] D. J. Gross and A. Neveu, Phys. Rev. D **10**, 3235 (1974).

[31] J.-L. Kneur and A. Neveu, Phys. Rev. D **85**, 014005 (2012).

[32] J.-L. Kneur and A. Neveu, Phys. Rev. D **88**, 074025 (2013).

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

[34] J.-L. Kneur and A. Neveu, Phys. Rev. D **92**, 074027 (2015).

[35] J.-L. Kneur and M. B. Pinto, Phys. Rev. D **92**, 116008 (2015).

[36] J.-L. Kneur and M. B. Pinto, Phys. Rev. Lett. **116**, 031601 (2016).

[37] G. N. Ferrari, J.-L. Kneur, M. B. Pinto, and R. O. Ramos, Phys. Rev. D **96**, 116009 (2017).

[38] P. Forgács, F. Niedermayer, and P. Weisz, Nucl. Phys. B**367**, 123 (1991).

[39] J. Gracey, Int. J. Mod. Phys. A **09**, 567 (1994); Phys. Lett. B **297**, 293 (1992).

[40] B. Kastening, Phys. Rev. D **54**, 3965 (1996); **57**, 3567 (1998).

[41] J.-L. Kneur, M. B. Pinto, and R. O. Ramos, Phys. Rev. D **74**, 125020 (2006).

[42] S. K. Gandhi, H. Jones, and M. B. Pinto, Nucl. Phys. B**359**, 429 (1991).

[43] I. T. Drummond, R. R. Horgan, P. V. Landshoff, and A. Rebhan, Nucl. Phys. B**524**, 579 (1998).

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

[45] J. Vermaseren, S. Larin, and T. van Ritbergen, Phys. Lett. B **405**, 327 (1997); M. Czakon, Nucl. Phys. B**710**, 485 (2005); K. Chetyrkin, Nucl. Phys. B**710**, 499 (2005).

[46] J. I. Kapusta and C. Gale, *Finite-Temperature Field Theory: Principles and Applications*, Cambridge Monographs on Mathematical Physics (Cambridge University Press, Cambridge, England, 2011).

[47] J.-L. Kneur and M. B. Pinto (to be published).