Renormalization group improved pressure for hot and dense quark matter

We apply the renormalization group optimized perturbation theory (RGOPT) to evaluate the quark contribution to the QCD pressure at finite temperatures and baryonic densities, at next-to-leading order (NLO). Our results are compared to NLO and state-of-the-art higher orders of standard perturbative QCD (pQCD) and hard thermal loop perturbation theory (HTLpt). The RGOPT resummation provides a nonperturbative approximation, exhibiting a drastically better remnant renormalization scale dependence than pQCD, thanks to built-in renormalization group invariance consistency. At NLO, upon simply adding to the RGOPT-resummed quark contributions the purely perturbative NLO glue contribution, our results show a remarkable agreement with ab initio lattice simulation data for temperatures $0.25 \lesssim T \lesssim 1 \, {\rm GeV}$, with a remnant scale dependence drastically reduced as compared to HTLpt.

We apply the renormalization group optimized perturbation theory (RGOPT) to evaluate the quark contribution to the QCD pressure at finite temperatures and baryonic densities, at nextto-leading order (NLO). Our results are compared to NLO and state-of-the-art higher orders of standard perturbative QCD (pQCD) and hard thermal loop perturbation theory (HTLpt). The RGOPT provides an all order resummed pressure in a well-defined approximation, exhibiting a drastically better remnant renormalization scale dependence than pQCD, thanks to built-in renormalization group invariance consistency. At NLO, upon simply adding to the RGOPT-resummed quark contributions the purely perturbative NLO glue contribution, our results show a remarkable agreement with ab initio lattice simulation data for temperatures 0.25 < ∼ T < ∼ 1 GeV, with a remnant scale dependence drastically reduced as compared to HTLpt.

I. INTRODUCTION
The complete prediction of the phase diagram describing strongly interacting matter transitions represents one of the major theoretical challenges in contemporary particle physics, despite the enormous progress achieved by lattice QCD (LQCD) numerical simulations. The main reason is that the well documented sign problem [1], which arises when finite chemical potential (µ) values are considered, prevents LQCD to be reliably applied at intermediate to high finite baryonic densities, while at low densities the problem may be circumvented, e.g., by performing a Taylor expansion around vanishing chemical potential results. In particular, within the latter regime LQCD has been very successful in predicting [2] that a crossover occurs at a pseudocritical temperature close to T pc ≈ 155 MeV when µ = 0. One alternative to describe the low temperature-high density domain is to employ effective quark theories [3], or the Nambu-Jona-Lasinio (NJL) model [4], evaluating physical quantities within some analytical more nonperturbative framework (e.g., mean field theory, MFT). This approach predicts that the (chiral) phase transition at low-T and finite µ is of the first kind [5] so that, as a byproduct, one should observe a critical end point (CP) signalled by a second order phase transition taking place at intermediate values of T and µ where the first order transition boundary terminates. This intriguing possibility is about to be tested in heavy-ion collisions experiments by decreasing the beam energy, √ s N N , so that the baryonic density increases. In view of these experiments it is an unfortunate situation that theoretical predictions using the full QCD machinery cannot be consistently carried out with the currently available nonperturbative techniques. As already emphasized LQCD is plagued by the sign problem while analytical tools such as the large-N approximation (which is related to MFT) may produce misleading results at criticality. More analytical alternatives to LQCD can partly address the deconfinement and/or chiral symmetry restoration at finite temperature and/or density: typically, some extensions of the NJL model [5,6] or Polyakov NJL (PNJL) [7,8], or other approaches incorporating more basic QCD dynamics in well-defined approximations, like the Dyson-Schwinger equations (see e.g. [9,10]), the functional Renormalization Group (see e.g. [11]), or other approaches [12]. On another side, standard thermal perturbation theory (PT) is unreliable at the relevant temperature and chemical potential ranges. Indeed, despite the asymptotic freedom (AF) property, its convergence can only be achieved at temperatures many orders of magnitude larger than the critical one. Even at intermediate temperatures, it is well-known that thermal PT is plagued by severe infrared divergences from bosonic zero modes, and has to be resummed to be more compatible with strong coupling regimes (for pedagogical reviews and lecture notes see, e.g., Refs. [13,14] and the very recent Ref. [15]). Yet, even the state-of-the-art, highest available order thermal PT [16,17], that incorporates a suitable resummation of infrared divergences, becomes more poorly accurate at moderate to low temperatures. A very successful alternative resummation method is to systematically expand from the start about a quasiparticle mass [18][19][20], that also more directly avoids infrared divergences apart from improving convergence issues. (One should keep in mind however that not all infrared divergences can be cured by such resummations: as is well-known the screening of static magnetic fields [21] is an additional nonperturbative phenomenon occuring in Yang-Mills theory, intrinsically limiting in practice the weak coupling expansion to maximal order α 3 S ln α S ). The expansion about a quasiparticle mass is actually close to analytical resummation approaches also used at zero-temperature, reminiscent of the traditional Hartree approximation and its variational generalizations, suitable to tackle the infrared divergence issues of massless theories. Basically one essentially deforms the original Lagrangian by a Gaussian mass term to be treated as an interaction, defining a modified perturbative expansion leading to a sequence of (variationally improved) approximations at successive orders.
The latter approaches appear under various names in the literature, such as optimized perturbation theory (OPT) [22][23][24] (as we dub it here), linear δ expansion (LDE) [25], variational perturbation theory (VPT) [26], or screened perturbation theory (SPT) [19,27] in the thermal context. Remark that adding a Gaussian term does not change the polynomial structure of the theory so that the process is compatible with the usual renormalization procedure. Already at NLO one usually goes beyond the simple Hartree approximation since the variational mass is dressed by incorporating different resummed topologies (exchange terms, vertex corrections, etc) order by order. Moreover, at leading order the OPT has the welcome property of exactly reproducing large-N results [28]. As discussed, e.g., in Ref. [29] this technique has been used to describe successfully a variety of physical situations, involving phase transitions in different models. On the other hand, for thermal theories, the SPT method has been generalized over the past two decades in order to be compatible with Yang-Mills theories. This generalization was made possible thanks to the hard thermal loop (HTL) gauge-invariant effective Lagrangian originally built by Braaten and Pisarski [18], consistently embedding HTL contributions, Landau damping and a screening gluon thermal mass term, with momentum-dependent self-energies and HTL-dressed vertices. The high temperature expansion based on the HTL Lagrangian, known as hard thermal loop perturbation theory (HTLpt) [20], has been employed in a series of applications up to NNLO (three-loops), to describe the QCD thermodynamics, considering both the glue [30] and quark [31][32][33] sectors at finite temperatures and baryonic densities. Given the intrinsic technical difficulties associated with the HTLpt evaluations, the NNLO state-of-the-art calculations performed typically in Refs. [32,33] represents a remarkable achievement. Unfortunately it is worth noting a serious remnant issue, also plaguing standard PT but not sensibly reduced in HTLpt: namely the sensitivity to the arbitrary renormalization scale is observed to substantially increase when higher orders are considered. More precisely, as compared to PT the NNLO HTLpt predictions in Refs. [32,33] are very close to the lattice results for temperatures down to T > ∼ 2 T pc for the commonly chosen "central" renormalization scale choice, M = 2π T 2 + µ 2 /π 2 . However, even a moderate scale variation of a factor 2 dramatically affects the pressure and related thermodynamical quantities by relative variations of order 1 or more. It has been argued [33] that resumming logarithms may help to improve the situation but, as explained in Refs. [34,35], it appears that the lack of renormalization group (RG) invariance is more basically rooted within the HTLpt approach. More recently an alternative combining the OPT framework with RG properties has been proposed: the renormalization group optimized perturbation theory (RGOPT) [34][35][36][37]. The main novelty is that it restores RG invariance at all stages of the calculation, in particular when fixing the arbitrary variational mass parameter. At vanishing temperatures it has been used in QCD up to high (three and four-loop) orders to estimate the basic coupling α s [37], predicting values very compatible with the world averages [38]. Also accurate values of the (vacuum) quark condensate were obtained at four-loop [39] and five-loop [40] orders. Concerning thermal theories the RGOPT has been applied to the simpler scalar φ 4 model [34,35] at NLO, as well as to the nonlinear sigma model (NLSM) [41]. In these thermal applications the RGOPT and PT/SPT predictions for the pressure have been compared, showing how the former approximation ameliorates the generic residual scale dependence of thermal perturbation theories at increasing perturbative orders. More recently we have evaluated the quark contribution to the QCD pressure at two-loop (NLO) at finite densities and vanishing temperatures, showing how the method improves over perturbative QCD (pQCD) [42]. In the present work we extend our approach to include the effects of a thermal bath. Note that applying the RGOPT readily to the glue contributions is beyond the present scope, due to some specific technical difficulties as briefly explained below (work in this direction is in progress [43]). Therefore in the present application the RGOPT resummation will be applied strictly only to the quark sector, while the gluons will be treated as in standard (NLO) PT. In the end both contributions will be combined in order to produce our complete final prediction for the NLO QCD pressure.
The paper is organized as follows. In the next section we briefly review our starting point, the perturbative expressions considered for the (massive) quark pressure at NLO, for which the basic RGOPT construction is recalled. In Sec. III the RGOPT is precisely defined for the quark pressure up to NLO (two-loop). Details of our two possible prescriptions at NLO are described in Sec. IV (that may be skipped by the reader only interested in the main results). Then Sec.V illustrates our main results for the pressure, both for the pure quark sector and for the full QCD one. We compare our results with the NLO and state-of-the-art higher orders of both PT and HTLpt, and also to lattice data for the complete QCD pressure. Sec. VI contains our conclusions and perspectives. Finally, three appendices specify some formulas and additional details used in our analysis.

II. QUARK CONTRIBUTION TO THE QCD PRESSURE
A. RG invariant perturbative pressure At two-loop order-g 1 , the contribution of (massive) quarks to the QCD perturbative pressure (the Feynman diagrams displayed in Fig. 1) can be obtained by combining the vacuum (T = µ = 0) results of Ref. [39] and the T, µ = 0 results of Refs. [44,45]. Considering the case of degenerate masses m u = m d = m s ≡ m, the renormalized pressure in the MS renormalization scheme, normalized per flavor, is in terms of the Fermi-Dirac distributions for anti-quarks (+ sign) and quarks (− sign), where µ represents the quark chemical potential, which relates to the baryonic chemical potential via µ B = 3µ.
In the present work we consider the case of symmetric quark matter and so do not distinguish the chemical potentials associated with the different flavors (µ s = µ u = µ d ≡ µ). The generalization to the case of chemical equilibrium needed to impose, e.g., β equilibrium should be straightforward. Also relevant for our purpose and comparisons is the wellknown resulting two-loop pressure expression for strictly massless quarks (that simplifies considerably since the J i integrals reduce to simple analytical expressions in this case, given for completeness in Appendix A): (2.8) withμ = µ/(2πT ). The Stefan-Boltzmann (SB) ideal gas limit reads Coming back to the massive quark case, we next define the standard homogenous RG operator, where our normalization convention for the QCD β-function and anomalous mass dimension γ m is where to two-loop order, 14) As is well known, the massive pressure (equivalently the massive free energy) is lacking perturbative RG-invariance: it can be easily checked that applying Eq.(2.10) to Eq.(2.1) leaves a remnant term starting at leading order. Now to turn Eq.(2.1) into a (perturbatively) RG invariant quantity, we proceed as in Refs. [34,35,39] (or closer to the present case, as in Ref. [42]), by subtracting a finite zero-point contribution, where the s i are determined at successive perturbative orders so that up to neglected higher order terms. Since our evaluations are being carried up to NLO, to restore perturbative RG invariance at this order it is sufficient to add the first two s 0 , s 1 coefficients that involve the coefficients of β(g), γ m (g) through Eq.(2.10). One finds explicitly [39,42] B. Implementing the RGOPT for the quark pressure The RGOPT construction [35,37,42] may be summarized as the following successive steps: • 1. First one restores perturbative RG invariance of the massive theory, following the procedure above described, giving at NLO P RGP T (m, g) in Eqs.(2.16), (2.18), (2.19). The subtraction contributions can be viewed as extending the Lagrangian by vacuum terms ∼ (m 4 /g)H(g), that could be equivalently introduced at the bare level [34,40], although the above prescription working directly with renormalized quantities is most convenient. At this stage, upon using commonly defined NLO running coupling and masses from Eqs.(2.13)-(2.15), it could easily be checked that Eq.(2.16) has a remnant scale dependence only of higher order ∼ m 4 g 2 by construction.
• 2. Next the RGOPT requires to variationally deform the Lagrangian, by rescaling the coupling (consistently for every standard interaction terms), and modifying quark mass terms, as (2.20) m being now an arbitrary mass. Eq.(2.20) is to be viewed as a convenient bookkeeping prescription actually performed on an already renormalized quantity P (m, g), thus obtained from standard perturbative expansion and appropriate renormalizations of the massive theory. Concerning the pressure it amounts to do those substitutions within P RGP T (m, g) of Eq.(2.16), thus incorporating the previously obtained vacuum subtractions 2 .
• 3. The resulting expression from Eq.(2.20) is expanded in powers of δ, at the same perturbative order considered, and next setting δ → 1 to recover the massless limit. At any finite order this leaves a residual m-dependence, that can be fixed by a stationarity criterion [22], the mass optimization prescription (MOP):  [37] the perturbative asymptotic freedom (AF) behavior for g → 0 at T = 0. This simple but compelling criterion often selects a unique solution, even at fiveloop order so far explored [40], in contrast with the related OPT/SPT approaches where using solely Eq.(2.21) generates an increasing number of possible solutions at increasing orders.
We recall that the quark sector at finite temperature has no zero modes, thus strictly there is no need to resum infrared divergences, and in standard thermal PT this sector is commonly treated purely perturbatively [16,17]. Accordingly since our present construction is being performed solely for the quark sector, it is not dealing with thermal infrared divergences (that anyhow occur only at NNLO from the gluon sector). It nevertheless resums well-defined RG-induced higher order contributions, leading to rather drastic differences with NLO pQCD, as will be illustrated below. At lowest nontrivial order δ 0 the resulting RGOPT pressure is given, keeping all terms of formally one-loop order, by Remark that the LO s 0 coefficient, Eq.(2.18), has produced the last term ∝ 1/b 0 in Eq.(2.24) after algebraic simplifications. There is a subtlety here: as Eq.(2.19) shows, s 1 involves two-loop RG coefficients and thus it is not mandatory to restore (perturbative) RG invariance at LO, that requires only s 0 = 0 as explained. Yet, since s 1 enters the pressure formally at O(1), it appears sensible to include it also within our one-loop RGOPT result Eq.(2.24), incorporating in this way higher order RG properties. (Actually the difference between the LO prescriptions with s 1 = 0 or taking more simply s 1 = 0 is not drastic). At the one-loop level the coupling runs according to the well-known expression . (2.25) Proceeding similarly at the next RGOPT order, the NLO pressure reads, after performing steps 1-4 above, The exact two-loop (2L) running coupling, analogue of the one-loop Eq.(2.25), is obtained by solving for g(M ) the implicit relation (see, e.g., Ref. [46]) for a given Λ MS value (this also defines the Λ MS basic scale in our normalization conventions). In the numerical illustrations below, we will use a value very close to the latest world average value [38], Λ MS = 335 MeV for N f = 3, that equivalently corresponds to α s (N f = 3, 1.5 GeV) 0.326. (NB the latter α s value precisely compares with the one taken in the literature for the NLO PT and HTLpt pressures [31]).

III. RG-OPTIMIZED RESUMMATION
A. One-loop RGOPT Before proceeding to our most relevant NLO results, derived basically from Eq.(2.26), it is useful to examine the probably more transparent RGOPT features at the lowest nontrivial (δ 0 ) LO. We recall that at this order the pressure already satisfies the massless RG Eq.(2.22) exactly, via the RG-driven exponent Eq.(2.23) of the variationally modified Lagrangian, Eq.(2.20). Consequently the arbitrary mass m may be fixed only by using the MOP Eq.(2.21). The latter acting on the LO pressure Eq.(2.24) can easily be cast into the form whose nontrivial solution gives an RG invariant dressed mass m(g, T, µ), since the combination 1/(b 0 g(M ))+ln m 2 /M 2 is trivially M -independent according to Eq.(2.25). (NB for more generality we keep s 1 unspecified at this stage, while for numerics below we will take s 1 = 0 as given by Eq. (2.19)). Once inserting m in Eq.(2.24) it produces a (one-loop) exactly RG invariant pressure, that takes the compact form: where it is understood that m is the nontrivial solution of Eq.(3.1) 4 . Some properties of the dressed mass m(g, T, µ) may be more transparent from considering the above expressions in the high temperature approximation (and µ = 0 to simplify), upon using well-known T m limits of the thermal integrals [14] J 1 , J 2 , given in Appendix A. This gives from Eq.(3.1) or, equivalently using Eq.(2.25) where we used 8π 2 s 1 = −53/84 for N f = 3. As seen in Eq.(3.3), for small coupling m admits a perturbative expansion having the expected form of a thermal screening mass. We stress however that m is unrelated to the perturbative Debye mass [44], which at one-loop order has the well-known expression (for µ = 0): In contrast m in Eq.(3.3) represents an intermediate variational quantity, whose meaning is merely once being inserted in P (m, g, T, µ) to define the (optimized) physical pressure at a given order. Remark that, upon embedding RG invariance properties via the subtraction terms in Eq.(2.16), leading to m in Eq.(3.1), the LO RGOPT pressure (3.2) involves nontrivial interaction terms. Indeed upon perturbatively reexpanding Eq.(3.2) using Eq.(3.3), it can be seen to resum arbitrary higher order contributions, although only those contributions induced by the specific leading order RG dependence 5 . Accordingly at LO and in the high-T approximation, using Eq.(2.25), Eq.(3.2) takes the simpler form Notice that the explicit dependence upon s 1 cancelled in P RGOP T 0 Eq.(3.2) upon using Eq.(3.1), but the solution m of Eq.(3.1) does depend on s 1 as specified. 5 In the simpler O(N ) φ 4 model, the analogous LO RGOPT [34] resums all large-N contributions, reproducing the exactly known large-N pressure [47], including nonanalytic terms ∼ λ 3p/2 , p ≥ 1, typical of a boson gas pressure.
normalized to the SB ideal quark gas P SB Eq.(2.9) (here for µ = 0). The fact that the higher order contributions may be absorbed essentially into a one-loop running coupling (for µ = 0 and high-T limits) is a peculiar LO feature of our construction: as we will see below at NLO the more involved RG-induced higher order corrections are not so simply incorporated. Another RGOPT feature is manifest in Eq.(3.6): at high-T the explicit M -dependence in Eq.(3.3) has been automatically traded for a dependence in g(∼ πT /Λ MS ), consequently from scale invariance, rather than being an extra convenient scale choice M ∼ πT to absorb ln(M/πT ) terms like in more standard (non-resummed) thermal perturbative expansions. The LO pressure Eq.(3.2) is however not expected to give a very realistic approximation of the complete higher order pressure, as it only relies on LO RG-invariance properties embedded within an essentially free gas pressure. The LO dressed mass m of Eq.(3.1) with exact T -dependence is illustrated as function of T in Fig. 3 (where it is also compared to NLO RGOPT dressed masses to be specified below). The corresponding pressure Eq.(3.2) is illustrated e.g. in Figs.5 and 6 for µ = 0 or in Figs. 9 and 10 for µ = 0, where it is also compared with NLO RGOPT and other NLO results. We will next proceed to the more realistic NLO RGOPT pressure: most of the above features will be maintained, except that the scale invariance can only be achieved approximately beyond LO, as we will examine.

B. Two-loop RGOPT
At NLO the RG Eq.(2.22) is no longer automatically satisfied by Eq.(2.26) with (2.23), and can be thus considered as an independent constraint. Following [35,41,42] we can in principle use either the MOP Eq.(2.21) or the RG Eq.(2.22), defining two possible alternative dressed mass m(g, T, µ): we will consider in the following both prescriptions, for completeness and comparison purposes. Accordingly the coupling g(M ) is simply determined from standard PT, i.e. with its running at (exact) two-loop order given by Eq.(2.27) and scale M chosen as a combination of πT and µ when both T, µ are non-zero. As shown above in Sec. III A, the LO RGOPT pressure exhibits one-loopexact scale invariance as a consequence of the simple form of Eq.(2.24) only depending on b 0 . While at NLO or beyond, P RGOP T (m(g), M, · · ·) inevitably has a remnant scale dependence: the basic reason is that the subtractions in Eq.(2.16) solely guarantee RG invariance up to remnant higher order terms, ∼ m 4 g 2 at NLO. This feature cannot be drastically reduced by the subsequent variational modification in Eq. (2.20). Concretely the exact NLO running coupling Eq.(2.27), that depends only on b i coefficients in Eqs.(2.13),(2.14), cannot perfectly cancel the explicit scale dependence of Eq.(2.26), that also involves anomalous mass dimension coefficients γ i in (2.15) reminiscent of the originally NLO massive theory. Nevertheless it is a nontrivial consequence of our above construction through steps 1)-4) in Sec.II B, preserving RG invariance at least at the same perturbative level as Eq. (2.17), that the remnant scale dependence remains formally of higher order ∼ m 4 g 2 . Accordingly this NLO RGOPT scale dependence, that we will exhibit by varying the scale by a factor 2 around central M ∼ 2πT (for µ = 0), is generically milder [34,35,41] than for standard PT and HTLpt. It is thus expected to remain moderate (and to further decrease at NNLO) even for relatively low temperature where the resulting dressed thermal mass is not necessarily perturbatively screened. Using the standard PT running coupling also more directly compares with the same common prescription in other related thermal resummations approaches, like HTLpt typically. But one should keep in mind that identifying the arbitrary renormalization scale M to be O(πT ) is strictly valid only at sufficiently high temperatures.
As mentioned above, Eq.(2.23) has the property to select a unique NLO solution matching the perturbative asymptotic freedom (AF) behavior for g → 0 at T = 0. As it happens however regarding the NLO quark pressure Eq.(2.26), imposing either Eq.(2.21) or Eq.(2.22) both fail to readily give a real dressed mass m(g, T, µ) for a substantial part of the physically relevant T, µ range. This is admittedly a technical burden of such methods, but the occurrence of complex variational solutions has no deeper physical meaning. Rather, it may be viewed to some extent as an accident of the specific MS scheme in which the original perturbative coefficients were calculated, given that nonreal solutions are often expected upon exactly solving nonlinear equations, like in the present case solving for m the NLO Eqs.(2.21) or (2.22). At the same time we wish to maintain these relations as exact as possible in order to capture RG resummation properties beyond PT. A crude escape could be simply to take the real part of the solutions, but that potentially loses some of the sought RG properties. The nonreal solution issue also occurred in the simpler T = µ = 0 case [37] as well as within the T = 0, µ = 0 cold quark matter application [42], where it was cured by performing a renormalization scheme change (RSC) [37]. The latter allows for the recovery of real solutions by modifying perturbative coefficients while keeping RG consistency by definition. Of course for such a solution to work the RSC should not be arbitrary, but fixed by a sensible prescription, and importantly such that it remains a moderate (i.e. perturbative) deviation from the original scheme. More specifically in [42] a relevant NLO RSC parameter B 2 was uniquely fixed by requiring collinearity of the RG and MOP curves in the {m, g} plane (that precisely expresses the recovering of real solutions). Technically this implies to nullify the determinant of partial derivatives of the RG and MOP equations, and to solve the latter together with, e.g., Eq.(2.21) for {B 2 , m(B 2 , g)}. While solving such a coupled system was easily manageable for the (entirely analytical) T = 0, µ = 0 NLO expressions in [42], it becomes numerically quite challenging for the rather involved T, µ = 0 NLO dependence from Eq. (2.26). Therefore in the present study, seeking as much as possible for simplicity, we will exploit the RSC arbitrariness quite similarly to recover real solutions, but via simpler alternative prescriptions precisely defined in next Sec. IV. The reader mainly interested in concrete results for the thermodynamical quantities may skip this section proceeding directly to Sec.V.

A. Simple RSC parametrization
Let us first specify for our later purposes the RSC to be used. Since one basically introduces a variational (quark) mass, the most natural and simplest RSC can be defined by modifying only the mass parameter: where a single B 2 coefficient parametrizes a perturbative NLO scheme change from the original MS-scheme 6 . As is well-known, for a perturbative series truncated at order g k (like in the present case the original order-g pressure Eq.(2.1)), different schemes differ formally by remnant terms of order O(g k+1 ), such that the difference between two schemes is expected to decrease at higher orders for sufficiently small coupling value. Note that we perform the perturbative RSC Eq. where for T, µ = 0, B mop and D mop take a relatively compact form: , (note that here x ≡ m 2 /T 2 ). In Eq.(4.4) we explicitly separated the T, µindependent part within D mop in the very first line to make its T, µ → 0 limit clear (remark also that D mop (T = 0) does not depend on m One easily recognizes that, for T → 0 the leading term for g → 0 in m 2 (−) has the correct AF behavior: ln m 2 M 2 (−) ∼ −1/(b 0 g), noting that b 0 = 9/(16π 2 ) (for N f = 3), which as recalled above is a compelling requirement of the RGOPT. In contrast the other (+) solution has a wrong sign and coefficient, thus drastically in contradiction with AF for g → 0. Therefore clearly only the above equation (4.2) with (−) is to be selected. It is further instructive to investigate the behavior of those two solutions for T = 0, taking for simplicity the high-T approximation (and µ = 0, see Eq.(A2)). After straightforward algebra one obtains, for the first few perturbative expansion terms: where we defined for short As seen the AF-compatible solution m(−) has a typical perturbative thermal screening mass behavior m ∼ √ g T , with a coefficient here mainly determined by RG properties (notice that the first order term is consistent with our LO above result, Eq.(3.3)). In contrast the non-AF-compatible Eq.(4.2) with (+) has m(+) solutions for T, µ = 0 having a coupling dependence that cannot be cast into the form of a perturbative expansion for small enough g. Moreover the corresponding real solutions generally give m/T 1, unless g is very large (see Appendix B). The latter features give further compelling reasons to rejecting this non-AF solution also for the T, µ = 0 case. Thus as anticipated the AF-compatibility criterion leads to a unique MOP solution.
The purely perturbative expansion Eq.(4.7) is however not expected to give a very good approximation for relatively low 8 T , and obviously not useful anyway for µ = 0. Before we proceed below with the more elaborate RSC prescription to solve exactly Eq.(4.2), it may be instructive to illustrate the results of using the simple perturbative solution Eq.(4.7), inserted in our NLO pressure Eq.(2.26), and with the resulting expression being truncated simply at first order in g (i.e. this is accordingly the NLO generalization of Eq.(3.6)). This is shown in Fig. 2, compared to the true standard NLO massless quark PT pressure Eq.(2.8). This result represents a good consistency check of our procedure: the two pressures are not strictly identical but very close, since after expressing the optimized mass m(g), the RGOPT is expected to approximate a massless theory. (Note that replacing m in Eq.(2.26) instead by, e.g., the standard thermal Debye mass Eq.(3.5), would give results more drastically departing from the massless PT pressure). Now more interestingly, the main purpose of the RGOPT is rather to provide higher order deviations from standard PT, induced by higher order RG-induced terms, as we will exhibit next.

C. NLO mass optimization prescription
Going back to the exact MOP Eq.(4.2), Eq.(4.4) involves the RSC parameter B 2 as induced from Eq.(4.1). In the original MS scheme, i.e. B 2 ≡ 0, D mop from Eq.(4.4) can take negative values for not particularly large couplings 9 . As anticipated above it therefore renders the (exact) m(g, T, µ) solution not always real, except in a rather limited range of physically interesting T and/or µ values. Remark however that since the (perturbatively leading) first term in Eq. results were obtained from modifying perturbative NLO expressions, one may simply expand perturbatively Eq.(4.2), obtaining therefore a real expression at arbitrary orders (as partially illustrated by the first few orders of such an expansion in Eq.(4.7)). But it is soon realized that this is a poor approximation of the actual exact expression, even for g slightly below the value at which D mop becomes negative. Accordingly such perturbative expansion would partly lose the sought RG properties, due to RG-consistent contributions being perturbatively truncated. Now with D mop not too far from being positive, a more efficient way to recover real solutions is from an appropriately chosen B 2 value such that D mop > 0.
Let us thus define precisely our prescription for the MOP Eq.(4.2): in a first stage we fix the arbitrary RSC parameter B 2 in Eq.(4.4) such that D mop > 0. Next, the resulting modified AF-matching Eq.(4.2) with (−) is solved exactly (numerically) for m(g, T, µ), recovering real solutions for practically most relevant g values. Note that simply requiring D mop ≥ 0 does not give a unique prescription, but it happens to be rather constrained: first D mop = 0 is excluded, as it would spoil the crucial AF-compatibility of Eq.(4.5), that at least requires the LO (first) term of Eq.(4.4). On the other hand if D mop > 0 would be too large, the AF-matching (−) Eq.(4.2) would take too negative values no longer giving a real solution (i.e., it cannot cross the x-axis). Since the problem comes from some negative terms within Eq.(4.4), a prescription that appears minimal is to fix B 2 such as to cancel solely the largest (in magnitude) T, µ-independent negative term within Eq.(4.4), −(47/6)g. Explicitly that gives: The latter B 2 prescription is very simple, and the resulting m M OP solution remains real for practically all physically relevant g(T, µ) values, while still including nontrivial higher order corrections induced from all remnant terms of Eq. (4.2). Other slightly different B 2 -fixing prescriptions are possible for T, µ = 0, but a notable property is that for different B 2 choices, that imply different exact m(B 2 ) solutions, the resulting physical pressure P (m(B 2 ), B 2 , · · ·) Eq.(2.26) happens to be largely insensitive to those unphysical B 2 parameter choices provided that m(B 2 ) remains real. This welcome feature is to be traced to the underlying RSC properties, together with the further perturbative screening from Eq. To close this subsection, we illustrate in Fig. 3 the resulting dressed thermal masses as function of the temperature, both at LO from Eq.(3.1), and NLO from MOP Eqs.(4.2), (4.9). As already mentioned their behavior is essentially that of screening thermal masses, except that those are determined from RG properties. We also compare in Fig. 3 with the similar dressed thermal mass as obtained from the alternative RG prescription, giving Eqs.(4.10) with (4.13), as we will specify in next subsection. Correspondingly Fig. 4 illustrates the relevant RSC deviation B 2 g 2 in Eq.(4.1) resulting from Eq.(4.9) as function of T . As an important crosscheck, it shows that the departure from the original MS-scheme remains quite moderate.

D. Alternative NLO RG prescription
Alternatively, the other very relevant prescription, as anticipated in Sec.III B, is to consider the RG Eq.(2.22) instead of the MOP Eq.(4.2) to determine the dressed mass m(g, T, µ). Once expressed for ln(m 2 /M 2 ) it takes a similar quadratic form as Eq.(4.2), conveniently normalized as where explicitly and Now, similarly to the previous MOP Eq.(4.2), for B 2 = 0 one obtains generally nonreal solutions since in D rg some contributions happen to be negative. In contrast with Eq.(4.2) however, the crucial AF-matching for the RG solution is already guaranteed solely from the first term in (4.11), up to higher order terms. These features strongly suggest the prescription fixing the arbitrary RSC parameter B 2 as simply to fully cancel D rg : Eq.(4.13) determines B 2 trivially using Eq.(4.12), leading to a single real AF-compatible solution m RG determined from the first two terms of Eq.(4.10), the latter being still an implicit equation in m for T, µ = 0 via J 2 entering Eq.(4.11). Eq.(4.13) may appear a rather peculiar choice, but there happen to be very few other choices to recover a real RG solution. We stress that for any (MOP or RG) prescriptions the resulting m(B 2 ) is an intermediate variational parameter without much physical meaning outside its use in the pressure. Here the resulting m RG (B 2 ) still involves arbitrary higher order contributions, as well as nontrivial T, µ dependence via B rg in Eq. (4.11). Similarly as for the MOP above prescription, we have checked that for other B 2 choices, as long as being moderately different from Eq.(4.13), our numerical RG results for T, µ = 0 are not strongly dependent upon those choices. The dressed exact thermal mass m RG resulting from Eqs.(4.10),(4.13) is illustrated as function of the temperature in Fig. 3, and compared with the previously discussed LO mass from Eq.(3.1) and m M OP from Eqs.(4.2), (4.9). As seen the dressed masses are numerically quite different, but such differences in the two alternative NLO variational masses are drastically reduced within the physical pressure as will be illustrated below. The corresponding RSC deviation B 2 g 2 obtained from Eq.(4.13) is illustrated in Fig. 4 as function of T , and compared to the similar MOP B 2 g 2 from Eq.(4.9). Notice that despite the visible discrepancies between the two expressions, they are numerically not drastically different and both behave smoothly, except at very low T < ∼ 0.5 GeV: as already mentioned above the important feature is that the induced departure from the original MS-scheme remains moderate.

V. RGOPT PRESSURE RESULTS AT NLO
To obtain the full benefit from the RGOPT, in particular the optimally reduced scale dependence, a price to pay as a result of the variational approach is to first solve exactly numerically for the dressed mass (either from Eq.(2.21) or alternatively Eq.(2.22)), prior to its use in the RGOPT pressure at NLO, Eq.(2.26). Such a procedure is moreover complicated by the onus of complex solutions, cured by the appropriate RSC as specified above in Sec. IV. But the relevant NLO expressions (4.2) or alternatively (4.10) are reasonably simple and the numerical procedure is straightforward. Before illustrating the resulting exact NLO RGOPT pressure, we start this section with another intermediate (more perturbative) prescription, to show the gradual improvement typically concerning the remnant scale dependence.

A. A simple perturbative approximation
The simplest we can do to recover real solutions without going through RSC considerations as elaborated on previously in Sec.IV, while capturing at the same time more accurate T, µ dependence, is to expand m from the MOP Eq.(4.2) perturbatively to NLO O(g 2 ), but keeping the exact thermal integrals in the resulting expression. This gives after simple algebra 10 therefore still to be solved numerically as an implicit function since J i ≡ J i ( m T , µ T ). The above expression readily gives a real solution, and allows to consider µ = 0 within the thermal integrals (and within the running coupling as well) while still keeping a relatively simple "perturbative-like" expression. Inserting the solution of Eq.(5.1) into the RGOPT NLO quark pressure Eq.(2.26) (keeping also exact thermal integrals consistently in the latter), gives the results illustrated for µ = 0 in Fig. 5, compared with the standard NLO PT pressure Eq.(2.8), and also with the NLO HTLpt (quark) pressure. (NB for a consistent comparison with the latter at this stage, we have extracted the sole quark contributions within the complete QCD NLO HTLpt pressure, which is not a trivial separation as in the case of NLO pQCD. How to do this precisely is explained in Appendix C). As seen in Fig. 5 the RGOPT pressure with the (MOP or RG) m(g 2 ) approximation has a more pronounced decrease, i.e. a departure from the ideal gas limit, than the standard NLO PT (pQCD) quark pressure and than LO RGOPT for moderate and low T values, that is mainly traced to the higher order g 2 contributions in Eq.(5.1) or Eq.(5.2). Actually, it is rather closer to the higher orders standard pQCD pressure, as will be illustrated below, partly due to Eq.(2.26) and the thermal functions J i being kept exact. (If perturbatively reexpanded, the resulting pressure gets back closer to the NLO pQCD result). This is in contrast with the NLO HTLpt pressure, that remains very close to the ideal gas limit except at very low T as seen in Fig. 5 11 . In Fig. 5 the RGOPT pressure also exhibits a better renormalization scale dependence as compared with NLO pQCD (at least for T > 1 GeV), although this is only a moderate improvement. Very similar results are obtained for µ = 0, that we omit to illustrate. We will see below that the more elaborate untruncated RGOPT pressure, accounting for higher orders in m(g), has a more drastically improved scale dependence, which is a main expected RGOPT feature.

MOP prescription
The resummation properties of the NLO RGOPT become more evident when one compares it with the standard perturbative one (pQCD) at the same NLO. We illustrate (first for µ = 0) the exact NLO RGOPT pressure P (m, g, T, µ) obtained from our first m M OP prescription, defined by solving Eqs.(4.2),(4.9) (as explained in details Subsec.IV C). In Fig. 6 the pressure is displayed as function of the temperature, compared with the LO RGOPT and the standard NLO pQCD Eq.(2.8), for the scale dependence πT ≤ M ≤ 4πT . The reduction of scale dependence stemming from the now exact (untruncated) NLO RGOPT appears substantial (about a factor ∼ 2 improvement for e.g. T ∼ 1 GeV). The HTLpt NLO (quark) pressure [31] is also shown in the same figure for comparison. We observe that the (NLO) quark HTLpt pressure has a small residual scale dependence for most T values (which is partly a consequence of limiting it to the quark only contribution), but does not depart very much from the ideal gas limit, in contrast with the RGOPT pressure. This latter feature is similar concerning the complete QCD NLO HTLpt [31]), while a more drastic departure from the ideal gas is only obtained at NNLO for HTLpt [33].

Alternative RG prescription
Similarly to Fig. 6, we illustrate in Fig. 7 the exact NLO RGOPT pressure as obtained from the alternative m RG prescription defined from solving Eqs.(4.10) and (4.13) (explained in details in Subsec.IV D). As is seen the RGOPT reduction of remnant scale dependence is even more substantial than for the previous m M OP prescription. The efficient reduction of remnant scale dependence with respect to standard NLO pQCD is also shown more quantitatively in Fig.  8, illustrating the maximal scale variations, ∆P/P ≡ (P (M = 4πT )/P (M = πT ) − 1, for the different approximations as indicated. Despite the numerically quite different MOP and RG dressed mass (see Fig 3), the resulting physical pressures are much closer for the two prescriptions, except at very low T values (i.e., very large coupling). This is a We now consider a nonzero chemical potential values. Since the MOP (4.2), (4.9) and RG (4.10), (4.13) prescriptions are defined quite generically they can be readily applied to the more general T, µ = 0 case. As a representative physical value we illustrate our results for µ B = 1.2 GeV. For the renormalization scale variation range we take as is common π T 2 + µ 2 /π 2 ≤ M ≤ 4π T 2 + µ 2 /π 2 within the exact NLO running coupling Eq. (2.27). This gives the results for the pressure as a function of temperature as shown in Fig. 9 and Fig. 10 for the MOP and RG prescriptions respectively. As is seen, for this rather sizable µ B value the qualitative picture is very similar to the µ B = 0 case above: namely the remnant scale dependence reduction from RGOPT is drastic as compared to pQCD, and sensible departures with respect to both pQCD and HTLpt are obtained from resummation effects at relatively low temperatures. These results appear to support the robustness of the RGOPT for a more reliable exploration of hot and dense matter.  Fig. 9 with alternative NLO RG prescription.

D. Including glue contribution: confrontation to lattice results
In principle a rather similar RGOPT treatment of the pure glue sector should be possible, building on the hard thermal loop (HTL) originally proposed in [18], with a gauge-invariant (non-local) effective Lagrangian properly describing Landau damping and screening with a gluon (thermal) "mass" term. However, this requires technically the evaluation of presently unknown and quite involved thermal integrals. More precisely, the RG-restoring subtraction analogous of Eq.(2.16) for nonzero gluon mass m D requires to calculate exact two-loop HTL m 4 g α S contributions, rather than expanded in m 2 D /T 2 up to order α 5/2 S , as calculated e.g. in [30,48]. Such a calculation involves up to five-dimensional complicated integrals, due to the highly nontrivial dressing of gluon propagators and vertices rooted in the HTL formalism. We leave such considerations for future work [43]. Therefore as above anticipated in the present work we treat the pure glue contribution most conservatively in a standard perturbative manner. At the same NLO, the standard perturbative pure glue contribution has the well-known expression [49] where the ideal gluon gas pressure is P g,SB = (8π 2 /45) T 4 . Thus we simply add the perturbative NLO contribution Eq.(5.3) (properly normalized) to our NLO RGOPT quark contributions Eq.(2.26), and for the numerical illustrations below we normalize our results to the full ideal pressure of quarks plus gluons: P SB → P q,SB + P g,SB

12
Following the progressive elaboration levels as in the previously shown quark pressure approximations, we first illustrate in Fig.11 the results of using the simple perturbatively re-expanded approximation for m, Eq.(5.1), for the quark contribution, but supplemented now by the NLO glue contribution, Eq.(5.3). The resulting RGOPT pressure is compared with both the (massless quark) state-of-the-art N 3 LO pQCD, which expression is taken from Ref. [16], and to available LQCD results from Ref. [50][51][52]. As is seen, adding the NLO PT glue contribution puts our results in the right ballpark of LQCD data, with clearly visible improvement as compared to pQCD, both for the central scale choice and resulting remnant scale uncertainty. (We also note that using instead the similar RG perturbative approximation Eq.(5.2) gives almost undistinguishable results from Fig. 11, illustrating the low order perturbative consistency of the two different MOP and RG prescriptions). Next in Figs.12 and 13, we illustrate similarly the results obtained upon adding the NLO PT glue contributions Eq.(5.3) to the NLO RGOPT quark pressure respectively for the (exact) MOP and RG prescriptions. These are compared with the state-of-the-art N 3 LO pQCD [16], and to LQCD results [50][51][52]. As seen the RGOPT results get closer to LQCD data, with a further reduced scale dependence, as compared to pQCD. In Fig.13 we compare in addition with both NLO [31] and the state-of-the-art NNLO HTLpt [33]. The corresponding HTLpt pressure expressions are worked out from Eqs.(51),(55),(56) in [31] at NLO, and from Eqs.(4.5),(4.6) in [33] at NNLO (we refer to Appendix C for more discussions on the HTLpt contributions). Notice also that these NNLO HTLpt and the O(g 3 ln g) pQCD [16]  Full NLO RGOPT (MOP prescription) plus NLO P P T g pressure as function of T (grey band) compared to (N 3 LO g 3 ln g) pQCD (light blue band), with scale dependence πT ≤ M ≤ 4πT , and to lattice data [50][51][52] at µB = 0.  [50][51][52] at µB = 0. results were obtained using a standard perturbative three-loop order running coupling. The pressure from the RG prescription gives the smallest residual scale uncertainties, and is in remarkable agreement with LQCD data in [50] for the central scale M = 2πT , for temperatures as low as T ∼ 0.25 GeV up to T = 1 GeV, the highest value considered in [50]. (More precisely let us mention that for the five available LQCD points in [50] with T > 0.3 GeV the central scale agreement is at the few permille level, and even slightly better when considering their estimated continuum data). It is also in good agreement with more recent LQCD data [51] at intermediate T . The RGOPT pressure is somewhat closer to LQCD results from [50] than the NNLO HTLpt pressure for 0.5 GeV < ∼ T < ∼ 1 GeV, while at higher T values HTLpt is nearer to the results of [52], and RGOPT shows more sizeable differences of order 5 − 7%. A concomitant feature however is the visible tension between low [50] and higher T [52] LQCD data in their common temperature range 13 .
Let us briefly mention that we have tried some variants of our prescriptions in order to check the stability of our results. First, the other RSC prescription to recover real solutions, mentioned above in Subsec. III B and used in Ref. [42], is to require the collinearity of the vectors tangent to the MOP and RG curves considered as functions of (m, g) (see Eq.(4.7) of Ref. [42]). In the present T = 0 case it is however numerically much more involved than our simpler prescriptions above (in particular to identify the AF-compatible solutions at moderate and low T values). Yet we could check that the resulting pressure is roughly similar to the one given by the MOP prescription in Figs. 6, 12. Next, we have also considered a variant of the RG prescription, by including the NNLO ∼ g m 4 s 2 subtraction term of Eq. (2.16), that is formally of NLO O(g) 14 . The s 2 expression [39,40] incorporates three-loop order RG coefficient dependence, thus for consistency we took a three-loop perturbative running coupling generalizing Eq. (2.27). We remark that the resulting pressure for this variant hardly shows visible differences with Figs. 13, reflecting a good stability, so that we omit to illustrate it.
Another physical quantity of interest is the trace anomaly (or equivalently interaction measure). The latter has the well-known expression (where the second and third equalities are of course valid only for µ = 0). As previously we add the pure glue NLO PT expression to our RGOPT quark contribution. The result is illustrated, for our best RG prescription, in Fig.14 where it is compared to LQCD data [50][51][52] only. A very good agreement with LQCD results of [50,51] is obtained for 0.3 GeV < ∼ T < ∼ 1 GeV, while there are more visible differences with the higher T results from [52]. Just for indication is also delineated the part of the remnant scale uncertainties originating solely from the RGOPT quark contributions (dashed lines) within the total uncertainties that also include the ones coming from the (standard) NLO PT glue contribution. As clearly seen however, similarly to pQCD and HTLpt, the NLO RGOPT does not describe correctly the peak region near the pseudocritical T c temperature as exhibited by lattice data. We speculate that a similar RGOPT resummation in the gluon sector should be a first necessary step to possibly better address the phase transition region, while at present we have treated this sector purely perturbatively as above explained. Therefore our present results are certainly not reliable in the region below T < ∼ 0.25 GeV. As a last fairly different illustration, we show the NLO RGOPT pressure as function of the quark chemical potential FIG. 14. NLO RGOPT (RG prescription) trace anomaly ∆ ≡ ε − 3P (including ∆ P T g ) (brown band) compared to lattice data [50][51][52]. The additional dashed lines illustrate the scale uncertainty originating solely from RGOPT quark contributions within the full scale uncertainty added by ∆ P T g (brown) band.
µ, for our two MOP and RG prescriptions respectively in Figs 15, 16, for a fixed relatively low temperature T = 0.3 GeV, thus definitely above the (µ = 0) pseudocritical temperature such that our present NLO approximation with above indicated limitations is presumably still reliable. As compared with HTLpt and pQCD, for both prescriptions the RGOPT pressure appears somewhat lower for high and moderate µ values, and exhibits a more regular behavior at low µ. The gain in remnant scale dependence appears once more drastic. We refrain to explore regions closer to the transition where our present construction becomes anyway unreliable. pressure (grey band), as function of the quark chemical potential µ for T = 0.3 GeV, compared to N 3 LO g 3 ln g pQCD (light blue band) and NNLO HTLpt (light red band), with scale dependence π(T 2 + µ 2 /π 2 ) 1/2 ≤ M ≤ 4π(T 2 + µ 2 /π 2 ) 1/2 .
To conclude this section it may be worth to recap the origin of the drastic differences between RGOPT and HTLpt, the latter being also basically a variational modification of the original QCD Lagrangian with mass terms, although based on the more elaborate HTL effective Lagrangian [18] (including among other features a thermal gluon mass parameter, m D ). There are essentially three important differences: • First, the perturbative RG-restoring subtraction terms, like in Eq.(2.16) typically, are missing in HTLpt. Accordingly the latter lacks perturbative RG-invariance formally by a leading order term of the massive theory pressure, O(m 4 ) ln(M/m). Now since for any (gluon or quark) thermal masses, m 2 ∼ #gT 2 , and HTLpt is also based on high temperature expansions, the latter uncancelled term is effectively only a three-loop order effect, thus largely screened and harmless at LO, and moderate even at NLO. In contrast this mismatch plainly resurfaces at NNLO HTLpt, presumably mainly explaining the large remnant scale dependence observed in Refs. [30,32,33].
• Second, the interpolating Lagrangian used in HTLpt is linear, namely with an exponent a = 1 in the HTL equivalent of Eq.(2.20), instead of our RG-determined Eq.(2.23). As we have shown [34] this generally spoils RG invariance even when the latter is fulfilled perturbatively by the original pressure.
• Finally, remark that upon choosing a variational mass prescription Eq.(2.21) in HTLpt (as was done e.g. in [31,32]), nonreal m may occur, similarly to what happens for RGOPT (although it happens rather at NNLO in HTLpt). In NNLO HTLpt applications this issue is avoided simply by replacing the gluon m D arbitrary mass by a perturbative thermal mass [30,33], and taking the quark mass m q = 0. However, enforcing perturbative masses is partly lacking the behavior beyond standard perturbation potentially provided by more variational prescriptions.

VI. CONCLUSIONS AND PERSPECTIVES
We have applied our RGOPT resummation approach at NLO at finite temperature and density for the QCD quark matter. As explained it generates more nonperturbative approximations with consistent RG properties already at LO (one-loop). Our NLO results have been compared to NLO and state-of-the-art N 3 LO pQCD predictions as well as to the state-of-the-art (NNLO) HTLpt results. Scale variations in the range πT ≤ M ≤ 4πT show that at NLO the method reduces scale uncertainties drastically as compared to pQCD. Since RG properties are consistently embedded within the RGOPT, we stress that generically the scale uncertainty bands observed at NLO should further shrink by considering the NNLO, O(g 2 ).
Our two possible 'MOP' or 'RG' prescriptions reflect the often non-uniqueness of variational approaches, although here their respective solution is unique from the compelling AF-matching requirement. Moreover the visible prescription difference for the resulting dressed mass (see Fig.3) is perturbatively consistent at low orders (Eqs.(5.1), (5.2)), and is substantially reduced within the resulting physical pressures. Using the RG Eq.(2.22) prescription, that more directly embeds consistent RG properties, not surprisingly gives the best remnant scale dependence at NLO (at it also happened in other considered models [34]). Once a specific RSC is adjusted to recover real solutions, the discrepancies between possibly different RSC prescriptions are formally perturbatively higher order terms. Nevertheless since we consider all expressions exactly rather than perturbatively truncated, numerically the RSC has a moderate net effect on the final pressure results. As we have illustrated, any perturbative reexpansion of the exact solutions somehow degrades the scale dependence.
Concerning the full QCD pressure, due to present technical limitations in applying the RGOPT plainly to the glue sector, in this work we have adopted a simple-minded approach, adding the purely perturbative NLO glue contributions to the pure quark sector resummed by RGOPT. We have confronted the resulting predictions for the QCD pressure with available LQCD results. For our best RG prescription the central scale M = 2πT results are in remarkable agreement with the LQCD results [50,51] for temperatures as low as T > ∼ 0.25 GeV, which lies within the nonperturbative regime, up to T = 1 GeV. However, similarly to pQCD and HTLpt, the NLO RGOPT construction explored in the present work is unable to describe the peak region near the pseudocritical T c temperature as exhibited by lattice data. Although our simple prescription appears to describe fairly well the moderate to high-T regimes T > ∼ 0.25 GeV ∼ 1.5 T pc , going beyond NLO one would not avoid to face the infrared divergences from gluon contributions, calling for appropriate resummations. The striking matching with LQCD results from Ref. [50] as seen in Fig.13 may be partly numerically accidental, but variants of our prescription, specifically the MOP pressure in Fig.12, still appears in very good agreement given our essentially NLO construction. Moreover the RG properties native to the RGOPT are not accidental in drastically reducing the scale dependence problem, particularly when comparing our NLO results to NNLO HTLpt. There are however some visible differences between our results and higher 1 GeV < ∼ T < ∼ 2 GeV LQCD data [52]. We remark that the LQCD pressure results in [50] and in Ref. [52] appear to be in tension in their common temperature range, while the trace anomaly shows more continuity 15 , a feature that may call for more investigations independently of our results. When comparing with 2 + 1 flavor LQCD as here illustrated, one may also keep in mind our presently not fully realistic approximation of N f = 3 degenerate flavors. As illustrated the RGOPT properties extend without much degradation to sizable chemical potential values and relatively low temperatures (see Figs. 9,10,and Figs. 15,16), that indicates the potential of our approach towards a more systematic exploration of hot and dense matter. Future applications may consider the inclusion of physical quark masses to generate a more realistic equation of state.
makes the AF solution identification obvious. Concerning the MOP Eq.(4.2), once B 2 is consistently determined by Eq.(4.9) such as to recover D mop > 0 in Eq.(4.4), one sees from the structure of (4. function of m/T and g(T /Λ M S ). Despite the nonlinear dependence in m/T , at the level of Eq.(4.2) both the AF and non-AF solutions happen to be unique in their respective existence range. This is illustrated in Fig. 17 (for µ = 0) for two representative low to moderate temperatures, respectively T = 0.5 and T = 1 GeV, and for the central scale choice M = 2πT . It is also clear that for any T the smallest solution is the AF one: Indeed for g(πT ≤ M ≤ 4πT ), − ln(m 2 /M 2 ) + B mop is a monotonically decreasing function of m for fixed T , and is > 0 (respectively < 0) below (respectively above) a given m 0 , such that necessarily m(AF) < m 0 < m(non-AF). The value of m 0 depends quite strongly on T (and M ): typically for the input corresponding to Fig. 17 with M = 2πT , one finds m 0 1.28(1.91) for T = 0.5(1.0) GeV respectively. (Notice also that in Fig. 17 the non-AF solution is unrealistically large with respect to T , that also makes it easy to unambiguously select the correct AF-matching solutions). At µ = 0, following the AF-matching m of Eq.(4.2) continuously from T = 0 to arbitrary T is in principle possible, although only for a fixed scale M (thus a fixed g(M )) unrelated to T , otherwise obviously at some small M ∼ πT one hits on M ∼ Λ M S where the perturbative coupling diverges. For sizable µ = 0 the latter problem if avoided if defining as conventional M ∼ π T 2 + µ 2 /π 2 (provided that one is not in the case of both T µ and small µ). Finally concerning the RG Eq.(4.10), both NLO solutions are already AF-matching, giving thus a unique solution upon using the prescription Eq.(4.13). Numerically the exact m RG solution of Eq.(4.10) is somewhat larger than m M OP for a given T , as illustrated in Fig. 3.

Appendix C: NLO and NNLO HTLpt expressions
For completeness we specify here how the NLO [31] and NNLO [33] HTLpt pressure expressions were precisely used when compared with other results. In particular for consistent comparison purposes in Figs. 5,6, and Figs. 9,10 we aim to pin down the HTLpt equivalent of the sole quark contributions, as shown up to NLO in Fig. 1, but with the quark and gluon propagators and quark-gluon vertex replaced with HTL-dressed ones consistently. More precisely from first comparing Eq.(51) of [31] to the pure glue NLO HTLpt pressure (given e.g. in Eq.(4.8) of second Ref. in [30]), it is not difficult to single out all terms originating solely from the pure quark vacuum energy. Next, from the resulting pressure we have rederived the (variationally determined) dressed thermal gluon m D and quark m q mass as in Eqs.(55),(56) in [31], that amounts to remove in these expressions the pure glue contributions (terms ∝ c A in Eq.(55) in [31]). At NNLO of HTLpt, a well-defined separation between pure quark and pure glue contributions appear ambiguous as these become more entangled. When comparing with our complete QCD RGOPT pressure e.g. in Fig.13 and subsequent figures, we obviously took the complete QCD NNLO HTLpt pressure, as given e.g. in Eqs.(4.5),(4.6) of