Renormalization group improved pressure for cold and dense QCD

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 non-trivial non-perturbative approximation which is completely renormalization group invariant. At the next-to-leading order the comparison between the RGOPT and the pQCD 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 $\bar{\rm MS}$ renormalization scale than those obtained with pQCD. 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.


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 non-linear and non-perturbative 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 Scam (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 end point. 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, α 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][13][14] apart from being rederived in a way compatible with the more modern 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: α 3 s ln 2 α 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 neutron stars [13,[17][18][19][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 MS arbitrary scale (M ) to be 2µ (where µ 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][22][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 µ ∼ 1 GeV, which corresponds to a baryon density ∼ 10 2 ρ 0 where ρ 0 ∼ 0.16 fm −3 represents the nuclear mass density. Such dependence indicates that the eventual non-perturbative 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πT is used at µ = 0 [22], while exhibiting a very large variation from M = πT to M = 4π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π T 2 + µ 2 /π.
In the present work we consider an alternative resummation method which incorporates RG properties to evaluate the QCD pressure at T = 0 and finite µ values at the two loop level. This technique, which provides non-perturbative 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 δ 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 = µ = 0 to evaluate the basic scale Λ MS [31,32] (equivalently the coupling α 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 λφ 4 theory [35,36], as well as to the non-linear 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 non-zero 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 = µ = 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][35][36][37]. The GN model is described by the Lagrangian density for a fermion field with N components given by [30] The theory described by Eq. (2.1) is invariant under the transformation ψ → γ 5 ψ 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 non-perturbative 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 2 GN = gπ/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 where δ is a book-keeping parameter interpolating between the free massive (δ = 0) and interacting massless (δ = 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][27][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 = 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 = µ = 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 = µ = 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 = 0 and/or µ = 0) After renormalizing in the MS-scheme (which at this lowest order amounts to simply a vacuum energy counterterm), one obtains where L m = ln(m/M ) and M is the arbitrary MS renormalization scale. Next consider the RG operator with the normalization conventions for the RG coefficients [39]: . 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: M dE/dM = −m 2 N/(2π) = 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 1 , that lead to the RG invariant (RGI) observable [32] Now requiring Eq. (2.8) to satisfy the RG equations fixes the s 0 coefficient to The procedure is easily generalized most conveniently when taking higher perturbative order contributions into account by considering a perturbative subtraction 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 φ 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-δ k , one aims to recover formally the massless limit, δ → 1. But the latter (re)expansion leaves a remnant dependence on the (arbitrary) mass m at any finite δ k order, since the expression was initially perturbative. Indeed, applying the RGOPT replacements, at lowest δ 0 order, to Eq. (2.8) gives 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: in consistency with the massless limit being seeked out. This uniquely fixes the exponent as a generic result also for other theories [32,[34][35][36][37]. Moreover, the same value of a is taken also when considering higher orders of the δ-expansion, keeping the simple interpolating form of Eq. (2.2), since this exponent is universal (renormalization scheme independent).
Thus substituting a = γ 0 /b 0 into Eq. (2.11) leads to the RGOPT lowest order result 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 Apart from the trivial result m = 0, one obtains which is clearly non-perturbative and explicitly RG invariant. Substituting 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-δ 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-δ (two loop level) and by applying the MOP criterion one fixes the mass to which is not RG invariant. As for the 1/N expansion, the first non trivial contribution appears at order-N 0 (the large-N limit, LN) whose gap equation yields the well known non-perturbative mass gap At this point a remarkable property of the RGOPT procedure over LN and standard OPT should be clear: it does produce a scale invariant non-perturbative 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 → ∞ 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, γ 0 = b 0 ≡ 1 at large-N 2 . The LN limit is also reproduced when considering in-medium effects: for example for the φ 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 "superdaisy" graphs, associated with plasmon infrared divergent contributions). Yet for more involved theories such as QCD, one does not expect the lowest δ 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 inmedium QCD case. Once having implemented the relevant RG subtraction coefficients in Eq.(2.10), one performs the RGOPT modification from Eq.(2.2) using the universal a value, Eq. (2.13), expanding this to δ k -order consistently with the original perturbative order considered, and taking the massless limit δ → 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 Λ MS at T = µ = 0 [32]. At non-vanishing temperatures the dressed mass also acquires a thermal dependence, but keeping its RG properties (see Refs. [35][36][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 (m) and a dressed running coupling (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 → 0 at T = 0, a simple criteria that most often selects a unique solution, even at the highest (fourloop) 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 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 non-perturbatively 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πα s ), in the limit of vanishing temperatures and finite baryonic densities with µ s = µ u = µ d ≡ µ, 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 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 section 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 (non-local) effective Lagrangian properly describing a gluonic (thermal) "mass" term in the HTLpt approximation [21][22][23]. However, in the present work, which deals only with the T = 0 and µ = 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 δ → 1 limit is taken, after the δ-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.
Feynman diagrams contributing to the perturbative quark pressure up to order-g.
To order-g the relevant contributions are displayed in Fig. 1. By combining the results of Ref. [34] for the vacuum (µ = 0) contributions with those of Ref. [12] for the in-medium part one obtains the renormalized result where p F = µ 2 − m 2 is the Fermi momentum, L m = ln(m/M ), d A = N 2 c − 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: 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 = µ = 0), with the appropriate QCD β and γ m RG functions. In our normalization conventions the QCD β(g) and γ m (g) functions read and where the coefficients are [45] b 0 = 1 and (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][36][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-δ 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 3 : in agreement with previous RGOPT applications to QCD [32,34]. The LO RGOPT pressure, per flavor, can then be written as Next, considering the NLO RGOPT (i.e. taking δ → 1 in the δ 1 -order expansion of Eqs. (3.1), (3.3)) and after some algebra the modified pressure (per flavor) reads: where P RGOP T 0,f (µ) is given by Eq. (3.13). Again, assuming µ > 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 5 L µ = ln[(µ+p F )/M ]. The LO and NLO RGOPT pressure are now ready to be optimized to generate non-perturbative approximation results as shown in the next section.

IV. OPTIMIZATION PROCEDURE AND NUMERICAL RESULTS
A. One-loop (δ 0 ) RGOPT Considering first the lowest (δ 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 RG = 0 exactly (at this order). Consequently the arbitrary mass m may be fixed only by using the MOP optimization equation: where Λ MS = M e − 1 2gb 0 is the one loop QCD Λ MS scale. Thus one obtains a non-perturbative mass, proportional to Λ 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 µ-dependent dressed mass m(µ) (while the reduced RG equation is still automatically satisfied at this order for µ = 0). At finite densities Eq. (4.1) is a little more involved due to the non-linear m-dependence from p F (m) = µ 2 − m 2 in Eq. (3.13). Yet, after simple algebra the formal solution may be cast into a compact form: (4.4) 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. 5 Those cancellations are however specific to the one-and two-loop level: at higher orders ln m and ln(µ + p F ) appear independently, due to more "nested" divergences in the bare calculation [13]. 6 Eq. 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 = 0 the term −1/2 in Eq. (4.4) above is simply replaced by −1/2 − 8π 2 s 1 = 11/84 (for n f = 3), which is the prescription used in the numerics below.
Eq. (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 ≡ g(M ) and the renormalization scale M only appear in the combination 1/(2b 0 g) + L µ 1/(2b 0 g) + ln(µ/M ) + · · ·, where the dots designate M -independent terms. Therefore, recalling that the (exact) one-loop running is defined as for a reference scale M 0 , it is immediate that Eq. For small coupling, the optimum mass m admits a perturbative expansion m 2 ∼ µ 2 (constant × g + O(g 2 )), which has the expected form of an (in-medium) screening mass. Nevertheless, we insist at this point that 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 (m, g, µ), that defines the (optimized) physical pressure at a given order. In fact, except for very weak coupling, the first order expansion of m 2 does not give a very good approximation of the exact m(µ): 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 m 2 (g(M )) < µ 2 even for large g(M ) (therefore consistent with the basic assumptions of the in-medium contributions). The numerical solution m at LO RGOPT from Eq. (4.3) as a function of µ is illustrated in Fig. 2, which among other features evidently confirms its exact scale invariance properties.

B. NLO two-loop (δ 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 non-linear 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 m(g), taking then for g more conservatively the purely perturbative two-loop running coupling. Unfortunately only non real 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 where the B i parametrize arbitrary scheme changes from the original MS-scheme 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 → 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 4 s 0 B 2 at two-loop order-g (and simply renaming m → 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 m solution the closest to the original 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): where f M OP and f 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 2 g 2 1.
As a first simple illustration, let us consider only the vacuum contribution at µ = 0. Applying Eq. (4.7), in conjunction with the MOP Eq.(4.1) and taking for g the two-loop perturbative Eq. We can now numerically optimize the NLO RGOPT pressure Eq. (3.14) including the in-medium contribution with µ = 0, adopting the RSC additional prescription to recover real m(µ) solutions for arbitrary µ 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 m(g), taking for g the purely perturbative running coupling, together with the RSC equation (4.7) to fix B 2 (µ) at NLO. (Clearly the optimal RSC parameter B 2 will now be a nontrivial function of the chemical potential µ, 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 µ ≤ M ≤ 4µ besides the "central" scale M = 2µ. The exact two-loop (2L) running coupling, analogue of the one-loop Eq.(4.5), is obtained by solving for g(M ) the relation for a given Λ MS value (this also defines the (two-loop order) Λ MS in our normalization conventions). Equivalently to giving a Λ MS value one can give a g(M 0 ) at some reference scale, M 0 . Here, we have chosen α s (M 0 = 1.5 GeV) = 0.326 to compare precisely with the values adopted in Ref. [16]: with (4.8) this corresponds to Λ MS 0.335 GeV, a value indeed very close to the present world average [33].  In contrast the NLO m(µ) 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 scaledependent terms of higher order O(g k+1 ), if the original perturbative expression is available at order g k . Fig. 3 illustrates the corresponding values of the RSC parameter combination B 2 (µ)g 2 (µ), thus quantifying the departure from MS-scheme. One can see that RSC remains reasonably perturbative, although the value of |B 2 (µ)| needed to recover real solutions are increasing rapidly for smaller µ values for the lower renormalization scale M = µ (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 equivalent massless free gas quantities P f g and ρ f g . These quantities, per flavor, are respectively and RGOPT O(g) FIG. 4. The normalized pressure as a function of the chemical potential at the central scale M = 2µ. 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).
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] P pQCD In Fig. 4 we show the normalized pressure predicted by the different order approximations at the central scale choice M = 2µ as adopted in Ref. [16]. The first thing to remark is that the RGOPT produces a non-trivial result already at order-δ 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 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 µ ≤ M ≤ 4µ, as in Ref. [13] where the scale dependence of pQCD results at orders g and g 2 have been analyzed 8 . The results are compared in Fig. 5, where the RGOPT appears to moderately improve the scale uncertainty, at least in the range µ > ∼ 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 ∆P/P f g ≡ (P (M = 4µ) − P (M = µ))/P f g as function of µ, 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 µ > ∼ 0.9 GeV (giving ∼ 25% improvement e.g. for µ 2 GeV), while the NLO pQCD scale dependence appears somewhat smaller in the lower µ range 0.5 < ∼ µ < ∼ 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 µ > ∼ 1.5 GeV. In contrast the RGOPT scale uncertainty is clearly better than the NNLO pQCD one in the full relevant µ range. We remark however that the smaller remnant dependence of NLO pQCD within the low-µ window (0.5 < ∼ µ < ∼ 0.9 GeV) is merely a side effect of the NLO pQCD pressure dropping towards zero at lower µ values than the two other approximations, as is clear from Fig. 5. Indeed not surprisingly all three approximations exhibit a rapidly growing scale dependence for µ values approaching the region where P (M µ) rapidly drops towards zero 9 . But Fig. 6 also shows that the maximal remnant dependence reached at the respective µ 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(µ)), none of the approximations is much reliable in the nonperturbative region where P/P f g (M µ) 1 due to large coupling (note, e.g., that µ < 0.8 already corresponds to α S (M = µ) > 0.5 using Eq.(4.8)). We thus conclude that, within the µ range where all the approximations are very reliable perturbatively, the NLO RGOPT remnant scale uncertainty is moderately but clearly improved in relative comparison to both NLO and NNLO pQCD (considering also that standard pQCD at T = 0, µ = 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 g m 4 s 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 = 0 at NLO scarcely changes our results (in contrast with the LO pressure where s 1 = 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: 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 ∼ µ) to obtain ρ(µ) in Eq.(4.12) involves a contribution ∝ ∂ M g(M ) on top of the explicit derivative ∂ µ 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 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 looses 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: 10 (4.13) Now, inserting this m expression into the NLO RGOPT pressure expression Eq. (3.14), with the running g → 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 µ ≤ M ≤ 4µ for the two expressions. One sees the quite remarkable agreement for the central scale choice M = 2µ (more precisely with less than ∼ 1.5% differences for any µ > 0.6 GeV), while the RGOPT scale "uncertainty" range is still slightly better even for this rather crude approximation. Concerning the scale dependence 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 µ 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 µ 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 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 non-perturbative approximations with consistent RG properties in a region of the QCD phase diagram which is currently unavailable to LQCD simulations. 10 In Eq.(4.13) the factors 9/7, 21/64, √ 43 are simply n f = 3 values of the specific combinations of RG coefficients b i , γ i appearing in this expression. The explicitly scale-dependent term ln µ/M only appears at next g 2 -order. Note also that we eliminated the other solution, with + √ · · ·, as it violates the necessary consistency m 2 ≤ µ 2 even for moderate g, and also does not fullfill the perturbative matching, ln µ/m ∼ 1/(2b 0 g), for µ m.
Our results have been compared to the state-of-the-art pQCD predictions that include a α 3 s ln 2 α s contribution. We have confirmed in this in-medium application the generic property that at lowest one-loop order this technique already captures non-trivial 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 µ 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µ the RGOPT results appear to be a very good approximation as they show a much better agreement with the perturbative higher orders O(α 3 s ln 2 α s ) than pQCD at the same order. Scale variations in the range M = µ − 4µ 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 = 0, µ = 0 compared with standard perturbation and HTLpt [35][36][37]. But this is merely due to the fact that standard pQCD at T = 0, µ = 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, eg, 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 specially when considering high temperature effects [47].