# Renormalization group improved pressure for hot and dense quark matter

Jean-Loïc Kneur<EMAIL_ADDRESS>Laboratoire Charles Coulomb
(L2C), UMR 5221 CNRS-Université Montpellier, 34095 Montpellier, France  Marcus
Benghi Pinto<EMAIL_ADDRESS>Departamento de Física, Universidade
Federal de Santa Catarina, 88040-900 Florianópolis, SC, Brazil  Tulio E.
Restrepo<EMAIL_ADDRESS>Departamento de Física, Universidade
Federal de Santa Catarina, 88040-900 Florianópolis, SC, Brazil  CFisUC -
Center for Physics of the University of Coimbra, Department of Physics,
Faculty of Sciences and Technology, University of Coimbra, 3004-516 Coimbra,
Portugal

###### Abstract

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.

##  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
($\mu$) 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}\approx 155\,{\rm MeV}$ when
$\mu=0$. One alternative to describe the low temperature-high density domain
is to employ effective quark theories [3], such as the Nambu–Jona-Lasinio
model [4], evaluating physical quantities within some analytical
nonperturbative framework (e.g., mean field theory, MFT). This approach
predicts that the (chiral) phase transition at low-$T$ and finite $\mu$ 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 $\mu$ 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, $\sqrt{s_{NN}}$, 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. On the other extreme 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, and
has to be resummed to be more compatible with strong coupling regimes (for
pedagogical reviews and lecture notes see, e.g., Refs. [6, 7] and the very
recent Ref. [8]). Yet, even the state-of-the-art, highest available order
thermal PT [9, 10], 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 around a quasiparticle mass [12, 11, 13], that also more
directly avoids infrared divergences apart from improving convergence issues.
This 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) [14, 15, 16] (as we dub it here), linear
$\delta$ expansion (LDE) [17], variational perturbation theory (VPT) [18], or
screened perturbation theory (SPT) [12, 19] 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 “nonperturbatively” dressed by
incorporating different 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 [20]. As discussed, e.g., in Ref. [21]
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 [11].
The high temperature expansion based on the HTL Lagrangian, known as hard
thermal loop perturbation theory (HTLpt) [13], has been employed in a series
of applications up to NNLO (three-loops), to describe the QCD thermodynamics,
considering both the glue [22] and quark [23, 24, 25] 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. [24, 25] 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.[24, 25] are very close to the lattice results for
temperatures down to $T\gtrsim 2\,T_{pc}$ for the commonly chosen “central”
renormalization scale choice, $M=2\pi\sqrt{T^{2}+\mu^{2}/\pi^{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 [25] that resumming logarithms may help to
improve the situation but, as explained in Refs.[27, 26], 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)[28, 29, 27, 26]. 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
$\alpha_{s}$ [29], predicting values very compatible with the world averages
[30]. Also accurate values of the (vacuum) quark condensate were obtained at
four-loop[31] and five-loop[32] orders. Concerning thermal theories the RGOPT
has been applied to the simpler scalar $\phi^{4}$ model [27, 26] at NLO, as
well as to the nonlinear sigma model (NLSM) [33]. 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)
[34]. 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 [35]).
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

###  II.1 RG invariant perturbative pressure

Figure 1: Feynman diagrams contributing to the perturbative quark pressure up
to NLO ${\cal O}(g)$.

At two-loop order-$g$ 111In all what follows we normalize for convenience the
(running) coupling in the $\overline{\rm MS}$-scheme as $g(M)\equiv
4\pi\alpha_{S}(M)$., 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=\mu=0$) results of Ref.[31] and the
$T,\mu\neq 0$ results of Refs.[36, 37]. Considering the case of degenerate
masses $m_{u}=m_{d}=m_{s}\equiv m$, the renormalized pressure in the
$\overline{\rm MS}$ renormalization scheme, normalized per flavor, is

$\displaystyle\frac{P^{PT}_{1}}{N_{f}\,N_{c}}$ $\displaystyle=$
$\displaystyle-\frac{m^{4}}{8\pi^{2}}\left(\frac{3}{4}-L_{m}\right)+2T^{4}J_{1}\left(\frac{m}{T},\frac{\mu}{T}\right)-3g\frac{m^{4}}{2\left(2\pi\right)^{4}}C_{F}\left(L_{m}^{2}-\frac{4}{3}L_{m}+\frac{3}{4}\right)$
(1) $\displaystyle-$ $\displaystyle
gC_{F}\left\\{\left[\frac{m^{2}}{4\pi^{2}}\left(2-3L_{m}\right)+\frac{T^{2}}{6}\right]T^{2}J_{2}\left(\frac{m}{T},\frac{\mu}{T}\right)+\frac{T^{4}}{2}J^{2}_{2}\left(\frac{m}{T},\frac{\mu}{T}\right)+m^{2}T^{2}J_{3}\left(\frac{m}{T},\frac{\mu}{T}\right)\right\\}\;,$

where $g\equiv g(M)$, $L_{m}\equiv\ln(m/M)$, $M$ is the arbitrary
renormalization scale, $C_{F}=(N_{c}^{2}-1)/(2N_{c})$, $N_{c}=3$, and
$N_{f}=3$. The in-medium and thermal effects are included in the
(dimensionless) single integrals:

$J_{1}\left(\frac{m}{T},\frac{\mu}{T}\right)=\int\frac{d^{3}{\bf\hat{p}}}{\left(2\pi\right)^{3}}\left\\{\ln\left[1+e^{-\left(E_{p}+\frac{\mu}{T}\right)}\right]+\ln\left[1+e^{-\left(E_{p}-\frac{\mu}{T}\right)}\right]\right\\}\;,$
(2)

with ${\bf\hat{p}}\equiv{\bf p}/T$,
$E_{p}=\sqrt{{\bf\hat{p}}^{2}+\frac{m^{2}}{T^{2}}}$,

$J_{2}\left(\frac{m}{T},\frac{\mu}{T}\right)=\int\frac{d^{3}{\bf\hat{p}}}{\left(2\pi\right)^{3}}\frac{1}{E_{p}}\left[f^{+}(E_{p})+f^{-}(E_{p})\right]\;,$
(3)

as well as in the double integral (after performing exactly the angular
integration over $p\cdot q/(|p||q|)$)

$J_{3}\left(\frac{m}{T},\frac{\mu}{T}\right)=\frac{1}{(2\pi)^{4}}\int^{\infty}_{0}\int^{\infty}_{0}\frac{d\hat{p}\,\hat{p}\,d\hat{q}\,\hat{q}}{E_{p}E_{q}}\left\\{\Sigma_{+}\ln\left[\frac{E_{p}E_{q}-\frac{m^{2}}{T^{2}}-\hat{p}\hat{q}}{E_{p}E_{q}-\frac{m^{2}}{T^{2}}+\hat{p}\hat{q}}\right]+\Sigma_{-}\ln\left[\frac{E_{p}E_{q}+\frac{m^{2}}{T^{2}}+\hat{p}\hat{q}}{E_{p}E_{q}+\frac{m^{2}}{T^{2}}-\hat{p}\hat{q}}\right]\right\\}\;,$
(4)

where

$\displaystyle\Sigma_{\pm}$ $\displaystyle=$ $\displaystyle
f^{+}\left(E_{p}\right)f^{\pm}\left(E_{q}\right)+f^{-}\left(E_{p}\right)f^{\mp}\left(E_{q}\right),$
(5)

in terms of the Fermi-Dirac distributions for anti-quarks ($+$ sign) and
quarks ($-$ sign),

$f^{\pm}(E_{p})=\frac{1}{1+e^{(E_{p}\pm\frac{\mu}{T})}}\;,$ (7)

where $\mu$ represents the quark chemical potential, which relates to the
baryonic chemical potential via $\mu_{B}=3\mu$.

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
($\mu_{s}=\mu_{u}=\mu_{d}\equiv\mu$). The generalization to the case of
chemical equilibrium needed to impose, e.g., $\beta$ equilibrium should be
straightforward. Also relevant for our purpose and comparisons is the well-
known 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):

$\frac{P^{PT}_{1}(m\to
0)}{P_{SB}(T,\mu)}=1-\frac{25g(M)}{42\pi^{2}}\,\left(\frac{1+\frac{72}{5}{\hat{\mu}}^{2}+\frac{144}{5}{\hat{\mu}}^{4}}{1+\frac{120}{7}{\hat{\mu}}^{2}+\frac{240}{7}{\hat{\mu}}^{4}}\right)$
(8)

with ${\hat{\mu}}=\mu/(2\pi T)$. The Stefan-Boltzmann (SB) ideal gas limit
reads

$P_{SB}(T,\mu)=T^{4}N_{f}N_{c}\left(\frac{7\pi^{2}}{180}\right)\left(1+\frac{120}{7}{\hat{\mu}}^{2}+\frac{240}{7}{\hat{\mu}}^{4}\right)\;.$
(9)

Coming back to the massive quark case, we next define the standard homogenous
RG operator,

$M\frac{d}{dM}=M\,\partial_{M}+\beta(g)\,\partial_{g}-m\,\gamma_{m}(g)\,\partial_{m}\;,$
(10)

where our normalization convention for the QCD $\beta$-function and anomalous
mass dimension $\gamma_{m}$ is

$\beta\left(g\right)=-2b_{0}g^{2}-2b_{1}g^{3}+\mathcal{O}\left(g^{4}\right)\;,$
(11)
$\gamma_{m}\left(g\right)=\gamma_{0}g+\gamma_{1}g^{2}+\mathcal{O}\left(g^{3}\right)\;,$
(12)

where to two-loop order,

$\displaystyle\left(4\pi\right)^{2}b_{0}$ $\displaystyle=$ $\displaystyle
11-\frac{2}{3}N_{f},$ (13) $\displaystyle\left(4\pi\right)^{4}b_{1}$
$\displaystyle=$ $\displaystyle 102-\frac{38}{3}N_{f},$ (14)
$\displaystyle\gamma_{0}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi^{2}}\;,\;\;\left(4\pi\right)^{4}\gamma_{1}=\frac{404}{3}-\frac{40}{9}N_{f}\;.$
(15)

As is well known, the massive pressure (equivalently the massive vacuum
energy) is lacking perturbative RG-invariance: it can be easily checked that
applying Eq.(10) to Eq.(1) leaves a remnant term of leading order. Now to turn
Eq.(1) into a (perturbatively) RG invariant quantity, we proceed as in
Refs.[31, 27, 26] (or closer to the present case, as in Ref.[34]), by
subtracting a finite zero-point contribution,

$\frac{P^{RGPT}}{N_{c}N_{f}}=\frac{P^{PT}}{N_{c}N_{f}}-\frac{m^{4}}{g}\sum_{k}s_{k}g^{k}\;,$
(16)

where the $s_{i}$ are determined at successive perturbative orders so that

$M\frac{d}{dM}\left(\frac{P^{RGPT}}{N_{c}N_{f}}\right)=0\,,$ (17)

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 $\beta(g)$, $\gamma_{m}(g)$ through Eq.(10). One finds
explicitly [31, 34]

$s_{0}=-\left[(4\pi)^{2}(b_{0}-2\gamma_{0})\right]^{-1},$ (18)
$s_{1}=-\frac{1}{4}\left[\frac{b_{1}-2\gamma_{1}}{4(b_{0}-2\gamma_{0})}-\frac{1}{12\pi^{2}}\right].$
(19)

###  II.2 Implementing the RGOPT for the quark pressure

The RGOPT requires to variationally deform the Lagrangian, by rescaling the
coupling (consistently for every standard interaction terms), while also
adding a modified Gaussian interpolating (mass) term. Explicitly, following
the prescription [29, 34]

${\cal L}_{QCD}^{RGOPT}={\cal L}_{QCD}|_{g\to\delta
g}-m(1-\delta)^{a}{\overline{\psi}}_{q}\psi_{q},$ (20)

where ${\cal L}_{QCD}$ is the standard massless QCD Lagrangian and $m$ is an
arbitrary quark mass at this stage. This is equivalent to perform the
substitution $m\to m(1-\delta)^{a}$ and $g\to\delta g$ in $P^{RGPT}$, Eq.
(16), reexpanded in powers of $\delta$, and finally setting $\delta\to 1$ to
recover the massless limit. At any finite order this leaves a residual
$m$-dependence, that is appropriately fixed by a stationarity criterion[14],
the mass optimization prescription (MOP):

$\frac{\partial{P}^{RGOPT}}{\partial m}\Bigr{|}_{\overline{m}}\equiv 0\;.$
(21)

The next step is to fix the arbitrary exponent $a$ that we introduced in
Eq.(20), by expanding to leading order-$\delta^{0}$ and requiring the
resulting pressure to satisfy, since Eq.(21) is used, the reduced (massless)
RG equation:

$\left[M\partial_{M}+\beta(g)\partial_{g}\right]P=0\;.$ (22)

This procedure yields

$a=\frac{\gamma_{0}}{2b_{0}}\,,$ (23)

that only depend on the universal (scheme-independent) LO RG coefficients, in
agreement with previous RGOPT applications to which we refer for demonstration
[29, 26, 34].

Accordingly at lowest nontrivial order the resulting RGOPT pressure is given,
keeping all terms of formally one-loop order, by

$\frac{P^{RGOPT}_{0}}{N_{f}N_{c}}=-\frac{2m^{4}}{\left(4\pi\right)^{2}}\left(\frac{3}{4}-L_{m}\right)-m^{4}\,s_{1}+2T^{4}J_{1}\left(\frac{m}{T},\frac{\mu}{T}\right)+\frac{m^{4}}{\left(4\pi\right)^{2}g\,b_{0}}\,.$
(24)

Remark that the LO $s_{0}$ coefficient, Eq.(18), has produced the last term
$\propto 1/b_{0}$ in Eq.(24) after algebraic simplifications. There is a
subtlety here: as Eq.(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}\neq 0$ as explained. Yet, since $s_{1}$ enters the
pressure formally at ${\cal O}(1)$, it appears sensible to include it also
within our one-loop RGOPT result Eq.(24), incorporating in this way a priori
more complete RG properties. (Actually, the difference between the LO
prescriptions with $s_{1}\neq 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

$g\left(M\right)=\frac{1}{2b_{0}\,\ln(M/\Lambda_{\overline{\rm MS}})}.$ (25)

Proceeding similarly at the next RGOPT order, the NLO pressure reads (after
setting $\delta=1$)

$\displaystyle\frac{P^{RGOPT}_{1}}{N_{f}N_{c}}$ $\displaystyle=$
$\displaystyle-\frac{m^{4}}{8\pi^{2}}\left(\frac{3}{4}-L_{m}\right)+2T^{4}\,J_{1}\left(\frac{m}{T},\frac{\mu}{T}\right)+\frac{m^{4}}{\left(2\pi\right)^{2}}\left(\frac{\gamma_{0}}{b_{0}}\right)\left(\frac{1}{2}-L_{m}\right)$
(26) $\displaystyle+$ $\displaystyle
m^{2}\left(\frac{\gamma_{0}}{b_{0}}\right)T^{2}\,J_{2}\left(\frac{m}{T},\frac{\mu}{T}\right)+\frac{m^{4}}{\left(4\pi\right)^{2}b_{0}}\left\\{\frac{1}{g}\left(1-\frac{\gamma_{0}}{b_{0}}\right)+\left[\left(b_{1}-2\gamma_{1}\right)\pi^{2}-\frac{\left(b_{0}-2\gamma_{0}\right)}{3}\right]\right\\}$
$\displaystyle-$ $\displaystyle
3gC_{F}\frac{m^{4}}{2\left(2\pi\right)^{4}}\left(L_{m}^{2}-\frac{4}{3}L_{m}+\frac{3}{4}\right)$
$\displaystyle-$ $\displaystyle
gC_{F}\left\\{\left[\frac{m^{2}}{4\pi^{2}}\left(2-3L_{m}\right)+\frac{T^{2}}{6}\right]T^{2}J_{2}\left(\frac{m}{T},\frac{\mu}{T}\right)+\frac{T^{4}}{2}J^{2}_{2}\left(\frac{m}{T},\frac{\mu}{T}\right)+m^{2}T^{2}\,J_{3}\left(\frac{m}{T},\frac{\mu}{T}\right)\right\\}.$

The exact two-loop (2L) running coupling, analogue of the one-loop Eq.(25), is
obtained by solving for $g(M)$ the implicit relation (see, e.g., Ref. [38])

$\ln\frac{M}{\Lambda_{\overline{\text{MS}}}}=\frac{1}{2b_{0}\,g}+\frac{b_{1}}{2b_{0}^{2}}\ln\left(\frac{b_{0}g}{1+\frac{b_{1}}{b_{0}}g}\right),$
(27)

for a given $\Lambda_{\overline{\text{MS}}}$ value (this also defines the
$\Lambda_{\overline{\text{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[30], $\Lambda_{\overline{\rm
MS}}=335\,{\rm MeV}$ for $N_{f}=3$, that equivalently corresponds to
$\alpha_{s}(N_{f}=3,1.5\,{\rm GeV})\simeq 0.326$. (NB the latter $\alpha_{s}$
value precisely compares with the one taken in the literature for the NLO PT
and HTLpt pressures[23]).

##  III RG-optimized resummation

###  III.1 One-loop RGOPT

Before proceeding to our most relevant NLO results, derived basically from
Eq.(26), it is useful to examine the probably more transparent RGOPT features
at the lowest nontrivial ($\delta^{0}$) LO. We recall that at this order the
pressure already satisfies the massless RG Eq.(22) exactly, via the RG-driven
exponent Eq.(23) of the variationally modified Lagrangian, Eq.(20).
Consequently the arbitrary mass $m$ may be fixed only by using the MOP
Eq.(21). The latter acting on the LO pressure Eq.(24) can easily be cast into
the form

$\frac{1}{b_{0}\,g}+\ln\frac{m^{2}}{M^{2}}-1-16\pi^{2}s_{1}-8\pi^{2}\frac{T^{2}}{m^{2}}J_{2}\left({\frac{m}{T}},{\frac{\mu}{T}}\right)=0,$
(28)

whose nontrivial solution gives a RG invariant dressed mass
$\overline{m}(g,T,\mu)$, since the combination $1/(b_{0}\,g(M))+\ln
m^{2}/M^{2}$ is trivially $M$-independent according to Eq.(25). (NB for more
generality we keep $s_{1}$ unspecified at this stage, while for numerics below
we will take $s_{1}\neq 0$ as given by Eq.(19)). Once inserting $\overline{m}$
in Eq.(24) it produces a (one-loop) exactly RG invariant pressure, that takes
the compact form:

$\frac{P^{RGOPT}_{0}}{N_{f}N_{c}}=2T^{4}\,J_{1}\left(\frac{\overline{m}}{T},\frac{\mu}{T}\right)+\frac{T^{2}}{2}{\overline{m}}^{2}J_{2}\left(\frac{\overline{m}}{T},\frac{\mu}{T}\right)-\frac{\overline{m}^{4}}{32\pi^{2}}\;,$
(29)

where it is understood that $\overline{m}$ is the nontrivial solution of
Eq.(28) 222Notice that the explicit dependence upon $s_{1}$ cancelled in
$P^{RGOPT}_{0}$ Eq.(29) upon using Eq.(28), but the solution $\overline{m}$ of
Eq.(28) does depend on $s_{1}$ as specified.. Some properties of the dressed
mass $\overline{m}(g,T,\mu)$ may be more transparent from considering the
above expressions in the high temperature approximation (and $\mu=0$ to
simplify), upon using well-known $T\gg m$ limits of the thermal integrals [7]
$J_{1},J_{2}$, given in Appendix A. This gives from Eq.(28)

${\overline{m}}^{2}(g,T,\mu=0)=T^{2}\frac{\pi^{2}}{3}\left[\frac{1}{2b_{0}g}-\ln\left(\frac{Me^{\gamma_{E}}}{\pi
T}\right)-8\pi^{2}s_{1}\right]^{-1}\;\simeq\frac{3}{8}g\,T^{2}+{\cal
O}(g^{2}),$ (30)

or, equivalently using Eq.(25)

${\overline{m}}^{2}(T,\mu=0)=T^{2}\frac{\pi^{2}}{3}\left[\ln\left(\frac{\pi
T}{e^{\gamma_{E}-\frac{53}{84}}\,\Lambda_{\overline{\rm
MS}}}\right)\right]^{-1}\;,$ (31)

where we used $8\pi^{2}s_{1}=-53/84$ for $N_{f}=3$. As seen in Eq.(30), for
small coupling $\overline{m}$ admits a perturbative expansion having the
expected form of a thermal screening mass. We stress however that
$\overline{m}$ is unrelated to the perturbative Debye mass [36], which at one-
loop order has the well-known expression (for $\mu=0$):

$m_{PT}^{2}=\frac{g}{6}\,T^{2}+{\cal O}(g^{2}).$ (32)

In contrast $\overline{m}$ in Eq.(30) represents an intermediate variational
quantity, whose meaning is merely once being inserted in
$P({\overline{m}},g,T,\mu)$ to define the (optimized) physical pressure at a
given order. Remark that, upon embedding RG invariance properties via the
subtraction terms in Eq.(16), leading to $\overline{m}$ in Eq.(28), the LO
RGOPT pressure (29) involves nontrivial interaction terms. Indeed upon
perturbatively reexpanding Eq.(29) using Eq.(30), it can be seen to resum
arbitrary higher order contributions, although only those contributions
induced by the specific leading order RG dependence333In the simpler $O(N)$
$\phi^{4}$ model, the analogous LO RGOPT[26] resums all large-$N$
contributions, reproducing the exactly known large-$N$ pressure[39], including
nonanalytic terms $\sim\lambda^{3p/2}$, $p\geq 1$, typical of a boson gas
pressure. . Accordingly at LO and in the high-$T$ approximation, using
Eq.(25), Eq.(29) takes the simpler form

$\frac{P^{RGOPT}_{0}}{P_{SB}}\simeq 1-\frac{5}{14}\left[\ln\left(\frac{\pi
T}{e^{\gamma_{E}-\frac{53}{84}}\,\Lambda_{\overline{\rm
MS}}}\right)\right]^{-1}\;,$ (33)

normalized to the SB ideal quark gas $P_{SB}$ Eq.(9) (here for $\mu=0$). The
fact that the higher order contributions may be absorbed essentially into a
one-loop running coupling (for $\mu=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.(33): at high-$T$ the explicit $M$-dependence in
Eq.(30) has been automatically traded for a dependence in $g(\sim\pi
T/\Lambda_{\overline{\rm MS}})$, consequently from scale invariance, rather
than being an extra convenient scale choice $M\sim\pi T$ to absorb $\ln(M/\pi
T)$ terms like in more standard (non-resummed) thermal perturbative
expansions.

The LO pressure Eq.(29) 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 $\overline{m}$ of Eq.(28) 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.(29) is illustrated e.g. in Figs.5 and 6 for $\mu=0$ or in Figs. 9 and 10
for $\mu\neq 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.

###  III.2 Two-loop RGOPT

At NLO the RG Eq.(22) is no longer automatically satisfied by Eq.(26) from
Eq.(23), and can be thus considered as an independent constraint. Following
[27, 33, 34] we can in principle use either the MOP Eq.(21) or the RG Eq.(22),
defining two possible alternative dressed mass $\overline{m}(g,T,\mu)$: 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.(27)
and scale $M$ chosen as a combination of $\pi T$ and $\mu$ when both $T,\mu$
are non-zero. A drawback however is that NLO partly spoils the sought RG
invariance properties: while at LO, as shown previously the running coupling
Eq.(25) perfectly matches the RGOPT pressure Eq.(24) such as to produce exact
(one-loop) scale invariance, at NLO the more involved running coupling Eq.(27)
has no reason to exactly match the scale dependence of Eq.(26). Accordingly a
remnant scale dependence inevitably arises at NLO, as we will exhibit by
varying the scale by a factor 2 around central $M\sim 2\pi T$ (for $\mu=0$).
Nevertheless, this RGOPT scale dependence is generically milder[27, 26, 33]
than the one produced by standard PT, as will be further illustrated in the
present analysis. This happens essentially since perturbative RG invariance is
maintained by construction for an arbitrary variational mass, such that the
RGOPT residual scale dependence is 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 ${\cal O}(\pi T)$ is strictly valid
only at sufficiently high temperatures.

A known unwelcome feature of any related OPT/SPT approaches is the generally
increasing number of possible solutions of e.g. Eq.(21) at increasing orders.
In contrast in our construction Eq.(23) further guarantees that the only
acceptable solutions are those matching[29] the perturbative asymptotic
freedom (AF) behavior for $g\to 0$ at $T=0$: a simple but compelling criteria
that often selects a unique solution, even at five-loop order so far explored
[29, 31]. As it happens however regarding the NLO quark pressure Eq.(26),
imposing either Eq.(21) or Eq.(22) both fail to readily give a real dressed
mass $\overline{m}(g,T,\mu)$ for a substantial part of the physically relevant
$T,\mu$ 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
$\overline{\rm 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.(21) or (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=\mu=0$ case [29] as well as within the
$T=0,\mu\neq 0$ cold quark matter application [34], where it was cured by
performing a renormalization scheme change (RSC)[29]. 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 [34] 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.(21) for $\\{B_{2},\overline{m}(B_{2},g)\\}$. While solving
such a coupled system was easily manageable for the (entirely analytical)
$T=0,\mu\neq 0$ NLO expressions in [34], it becomes numerically quite
challenging for the rather involved $T,\mu\neq 0$ NLO dependence from Eq.(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.

##  IV NLO prescriptions

###  IV.1 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:

$m\to m^{\prime}(1+B_{2}g^{2})\,,$ (34)

where a single $B_{2}$ coefficient parametrizes a perturbative NLO scheme
change from the original ${\overline{\text{MS}}}$-scheme 444Eq.(34) has also
the welcome property that it does not affect the definition of the reference
QCD scale $\Lambda_{\overline{\text{MS}}}$, in contrast with a similar
perturbative modification acting on the coupling (see Ref.[29] for details)..
As is well-known, for a perturbative series truncated at order $g^{k}$ (like
in the present case the original order-$g$ pressure Eq.(1)), different schemes
differ formally by remnant terms of order ${\cal 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.(34) consistently on the original PT expression (1) prior to its
(nonperturbative) modification induced from Eq.(20) with the subsequent
$\delta$-expansion. The net RSC modification to the pressure is to add an
extra term, $-4g\,(m^{\prime})^{4}s_{0}B_{2}$, entering thus the resulting
exact NLO Eq.(21) or Eq.(22). Thus Eq.(34) modifies the latter equations
purposefully, now considering those equations as constraints for the arbitrary
mass $m^{\prime}$, after the (nonperturbative) modifications from Eq.(20)555To
avoid excessive notation proliferation, in what follows once having performed
the replacement implied by Eq. (34) we simply rename $m^{\prime}\to m$ the
variational mass to be determined from Eq.(21) or Eq.(22).. Accordingly
$B_{2}$ may be considered as an extra variational parameter, quite similarly
to $m$, thus to be fixed by a definite prescription as will be specified
below.

###  IV.2 AF-compatible NLO dressed mass solutions

To identify some relevant properties of the sought dressed
$\overline{m}(g,T,\mu)$ solutions we consider first the MOP Eq.(21) more
explicitly at NLO, thus applied to Eq.(26). It is convenient to formally solve
it in a first stage for $\ln[m^{2}/M^{2}]$, as that would give simply an exact
quadratic equation at $T=\mu=0$. Accordingly the two equations (that are
implicit in $m$ for $T,\mu\neq 0$) can be conveniently written, after
straightforward algebra, as

$-\ln\frac{m^{2}}{M^{2}}+B_{mop}\mp\frac{2\pi}{3g}\sqrt{D_{mop}}=0\;,$ (35)

where for $T,\mu\neq 0$, $B_{mop}$ and $D_{mop}$ take a relatively compact
form:

$B_{mop}=-\frac{7\pi^{2}}{9g}+\frac{5}{6}+4\pi^{2}\left(J_{2}^{\prime}+\frac{T^{2}}{m^{2}}J_{2}\right),$
(36) $\displaystyle D_{mop}$ $\displaystyle=$ $\displaystyle
9\frac{\pi^{2}}{4}-\frac{47}{6}g-g^{2}\left(\frac{35}{16\pi^{2}}+288\frac{\pi^{2}}{7}B_{2}\right)$
(37)
$\displaystyle+36\pi^{2}g^{2}\left(J_{2}^{\prime\,2}+\frac{T^{4}}{m^{4}}J_{2}^{2}\right)+9g(g-2\pi^{2})\left(\frac{T^{2}}{m^{2}}J_{2}-J_{2}^{\prime}\right)$
$\displaystyle+8\pi^{2}g^{2}\frac{T^{2}}{m^{2}}\left(3J_{2}-1\right)J_{2}^{\prime}-48\pi^{2}g^{2}\left(J_{3}^{\prime}+\frac{T^{2}}{m^{2}}J_{3}\right)\;,$

where $J_{i}\equiv J_{i}(m^{2}/T^{2},\mu/T)$,
$J_{i}^{\prime}\equiv\partial_{x}J_{i}(x)$, (note that here $x\equiv
m^{2}/T^{2}$). In Eq.(37) we explicitly separated the $T,\mu$-independent part
within $D_{mop}$ in the very first line to make its $T,\mu\to 0$ limit clear
(remark also that $D_{mop}(T=0)$ does not depend on $m$).
One first property of Eq.(35) is exhibited from expanding it perturbatively to
the first few terms. That gives

$\ln\frac{\overline{m}^{2}}{M^{2}}(-)\simeq-\frac{16\pi^{2}}{9g}+\frac{139}{54}+8\pi^{2}\frac{T^{2}}{\overline{m}^{2}}J_{2}+{\cal
O}(g),$ (38)

and

$\ln\frac{\overline{m}^{2}}{M^{2}}(+)\simeq\frac{2\pi^{2}}{9g}-\frac{49}{54}+8\pi^{2}J_{2}^{\prime}+{\cal
O}(g).$ (39)

One easily recognizes that, for $T\to 0$ the leading term for $g\to 0$ in
$\overline{m}^{2}(-)$ has the correct AF behavior:
$\ln\frac{\overline{m}^{2}}{M^{2}}(-)\sim-1/(b_{0}g)$, noting that
$b_{0}=9/(16\pi^{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\to 0$.
Therefore clearly only the above equation (35) with $(-)$ is to be selected.
It is further instructive to investigate the behavior of those two solutions
for $T\neq 0$, taking for simplicity the high-$T$ approximation (and $\mu=0$,
see Eq.(52)). After straightforward algebra one obtains, for the first few
perturbative expansion terms:

$\frac{\overline{m}^{2}_{(-)}}{T^{2}}=\frac{3}{8}g\left[1-\frac{3}{8\pi^{2}}g\left(3L_{T}+\frac{85}{36}\right)\right]^{-1}+g^{2}\,\left(\frac{67}{288\pi^{2}}+6J_{3}(0,0)\right)+{\cal
O}(g^{3}),$ (40)

where we defined for short

$L_{T}\equiv\ln\left(\frac{Me^{\gamma_{E}}}{\pi T}\right).$ (41)

As seen the AF-compatible solution $\overline{m}(-)$ has a typical
perturbative thermal screening mass behavior $m\sim\sqrt{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.(30)). In contrast the
non-AF-compatible Eq.(35) with $(+)$ has $\overline{m}(+)$ solutions for
$T,\mu\neq 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\gg 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,\mu\neq 0$ case. Thus as anticipated the
AF-compatibility criteria leads to a unique MOP solution.

The purely perturbative expansion Eq.(40) is however not expected to give a
very good approximation for relatively low 666In particular the exact NLO
$\overline{m}$ as obtained below can be such that $\overline{m}/T>1$ at
sufficiently low $T$ (see Fig. 3), somewhat invalidating the high-$T$
approximation. $T$, and obviously not useful anyway for $\mu\neq 0$. Yet,
before to proceed below with the more elaborate RSC prescription to solve
exactly Eq.(35), it may be instructive to illustrate the results of using the
simple perturbative solution Eq.(40), inserted in our NLO pressure Eq.(26),
and with the resulting expression being truncated simply at first order in $g$
(i.e. this is accordingly the NLO generalization of Eq.(33)). This is shown in
Fig. 2, compared to the true standard NLO massless quark PT pressure Eq.(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 $\overline{m}(g)$, the RGOPT is expected to approximate a
massless theory. (Note that replacing $m$ in Eq.(26) instead by, e.g., the
standard thermal Debye mass Eq.(32), 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.

Figure 2: Perturbatively re-expanded NLO RGOPT pressure $P(T,\mu=0)$ (red
band) compared with standard perturbative NLO pressure Eq.(8) (blue band),
with scale dependence $\pi T\leq M\leq 4\pi T$.

###  IV.3 NLO mass optimization prescription

Going back to the exact MOP Eq.(35), Eq.(37) involves the RSC parameter
$B_{2}$ as induced from Eq.(34). In the original $\overline{\rm MS}$ scheme,
i.e. $B_{2}\equiv 0$, $D_{mop}$ from Eq.(37) can take negative values for not
particularly large couplings777 For example for $T=\mu=0$ where only the first
three terms of Eq.(37) are nonvanishing, $D_{mop}\leq 0$ for rather moderate
$g\geq 2.64$, i.e $\alpha_{S}\geq 0.21$.. As anticipated above it therefore
renders the (exact) $\overline{m}(g,T,\mu)$ solution not always real, except
in a rather limited range of physically interesting $T$ and/or $\mu$ values.
Remark however that since the (perturbatively leading) first term in Eq.(37)
is positive, this lost of real $\overline{m}$ solutions arises solely when
considering the exact Eq.(35): now since all our results were obtained from
modifying perturbative NLO expressions, one may simply (re)expand Eq.(35) in
perturbation, that gives a real expression at arbitrary orders (as partially
illustrated by the first few orders of such an expansion in Eq.(40)). 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 it would partly lose the sought good RG properties, due
to RG-induced 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.(35): in a first
stage we fix the arbitrary RSC parameter $B_{2}$ in Eq.(37) such that
$D_{mop}>0$. Next, the resulting modified AF-matching Eq.(35) with $(-)$ is
solved exactly (numerically) for $\overline{m}(g,T,\mu)$, recovering real
solutions for practically most relevant $g$ values. Note that simply requiring
$D_{mop}\geq 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.(38), that at least requires the LO (first)
term of Eq.(37). On the other hand if $D_{mop}>0$ would be too large, the AF-
matching $(-)$ Eq.(35) 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.(37), a prescription that appears minimal is to
fix $B_{2}$ such as to cancel solely the largest (in magnitude)
$T,\mu$-independent negative term within Eq.(37), $-(47/6)g$. Explicitly that
gives:

$B_{2}=-\frac{329}{1728\pi^{2}\,g}\;.$ (42)

The latter $B_{2}$ prescription is very simple, and the resulting
$\overline{m}_{MOP}$ solution remains real for practically all physically
relevant $g(T,\mu)$ values, while still including nontrivial higher order
corrections induced from all remnant terms of Eq. (35). Other slightly
different $B_{2}$-fixing prescriptions are possible for $T,\mu\neq 0$, but a
notable property is that for different $B_{2}$ choices, that imply different
exact $\overline{m}(B_{2})$ solutions, the resulting physical pressure
$P(\overline{m}(B_{2}),B_{2},\cdots)$ Eq.(26) happens to be largely
insensitive to those unphysical $B_{2}$ parameter choices provided that
$\overline{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.(40): ${\overline{m}^{2}}\sim(3/8)gT^{2}+{\cal O}(g^{2})$:
as easily checked, $B_{2}$ from Eq.(34) only appears at higher order ${\cal
O}(g^{3})$ both in the perturbatively expanded $\overline{m}$ Eq.(40) and
corresponding re-expanded pressure from Eq.(26). In other words once $B_{2}$
is adjusted to recover a real $\overline{m}$ solution of Eq.(35), the
discrepancies between possibly different $B_{2}$ prescriptions are somewhat
hidden within perturbatively higher order terms.
To close this subsection, we illustrate in Fig. 3 the resulting dressed
thermal masses as function of the temperature, both at LO from Eq.(28), and
NLO from MOP Eqs.(35), (42). 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, to be specified in next
subsection.

Figure 3: Exact LO RGOPT thermal mass (dot-dashed) compared with exact MOP and
RG NLO thermal mass for $\pi\leq M\leq 4\pi T$ at $\mu_{B}=0$.

Correspondingly Fig. 4 illustrates the relevant RSC deviation $B_{2}g^{2}$ in
Eq.(34) resulting from Eq.(42) as function of $T$. As an important crosscheck,
it shows that the departure from the original $\overline{\rm MS}$-scheme
remains quite moderate.

Figure 4: RSC parameter $B_{2}g^{2}(M)$ for the MOP and RG prescriptions for
$\pi T\leq M\leq 4\pi T$ at $\mu_{B}=0$.

###  IV.4 Alternative NLO RG prescription

Alternatively, the other very relevant prescription, as anticipated in
Sec.III.2, is to consider the RG Eq.(22) instead of the MOP Eq.(35) to
determine the dressed mass $\overline{m}(g,T,\mu)$. Once expressed for
$\ln(m^{2}/M^{2})$ it takes a similar quadratic form as Eq.(35), conveniently
normalized as

$-\ln\frac{m^{2}}{M^{2}}+B_{rg}\mp\frac{8\pi^{2}}{g}\sqrt{\frac{2}{3}D_{rg}}=0\;,$
(43)

where explicitly

$B_{rg}=-\frac{1}{b_{0}\,g}+\frac{172}{81}-\frac{64}{81}\left(\frac{4g}{9\pi^{2}}\right)\,\frac{1}{1+\frac{4g}{9\pi^{2}}}+8\pi^{2}\frac{T^{2}}{m^{2}}J_{2},$
(44)

and

$D_{rg}=-\left(\frac{3}{7}B_{2}+\frac{11}{384\pi^{4}}\right)g^{2}-\frac{g}{27}\frac{(4g+81\pi^{2})}{(4g+9\pi^{2})^{2}}+g^{2}\,\frac{T^{4}}{m^{4}}J_{2}\left(J_{2}-\frac{1}{6}\right)-g^{2}\,\frac{T^{2}}{m^{2}}J_{3}.$
(45)

Now, similarly to the previous MOP Eq.(35), for $B_{2}=0$ one obtains
generally nonreal solutions since in $D_{rg}$ some contributions happen to be
negative. In contrast with Eq.(35) however, the crucial AF-matching for the RG
solution is already guaranteed solely from the first term in (44), 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}$:

$D_{rg}(B_{2})\equiv 0.$ (46)

Eq.(46) determines $B_{2}$ trivially using Eq.(45), leading to a single real
AF-compatible solution $\overline{m}_{RG}$ determined from the first two terms
of Eq.(43), the latter being still an implicit equation in $m$ for $T,\mu\neq
0$ via $J_{2}$ entering Eq.(44). Eq.(46) 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
$\overline{m}(B_{2})$ is an intermediate variational parameter without much
physical meaning outside its use in the pressure. Here the resulting
$\overline{m}_{RG}(B_{2})$ still involves arbitrary higher order
contributions, as well as nontrivial $T,\mu$ dependence via $B_{rg}$ in
Eq.(44). Similarly as for the MOP above prescription, we have checked that for
other $B_{2}$ choices, as long as being moderately different from Eq.(46), our
numerical RG results for $T,\mu\neq 0$ are not strongly dependent upon those
choices.

The dressed exact thermal mass $\overline{m}_{RG}$ resulting from
Eqs.(43),(46) is illustrated as function of the temperature in Fig. 3, and
compared with the previously discussed LO mass from Eq.(28) and
$\overline{m}_{MOP}$ from Eqs.(35), (42). 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.(46) is illustrated in Fig. 4 as function of $T$, and
compared to the similar MOP $B_{2}g^{2}$ from Eq.(42). Note that despite the
visible discrepancies between the two expressions, they are numerically not
drastically different and both behave smoothly, except at very low $T\lesssim
0.5$ GeV: as already mentioned above the important feature is that the induced
departure from the original $\overline{\rm 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.(21) or
alternatively Eq.(22)), prior to its use in the RGOPT pressure at NLO,
Eq.(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 (35) or alternatively (43) 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.

###  V.1 A simple perturbative approximation

The simplest we can do to recover real solutions without going through RSC
considerations as elaborated previously in Sec.IV, while capturing at the same
time more accurate $T,\mu$ dependence, is to expand $\overline{m}$ from the
MOP Eq.(35) perturbatively to NLO ${\cal O}(g^{2})$, but keeping the exact
thermal integrals in the resulting expression. This gives after simple algebra
888As an algebraic subtlety, one should first expand perturbatively the (AF-
matching) Eq.(35) with $(-)$ before to solve it formally for $m^{2}/T^{2}$,
otherwise one loses the latter AF-matching properties.

$\frac{m^{2}_{MOP}}{T^{2}}=\frac{9}{2}gJ_{2}+g^{2}\left[\frac{17}{9}J^{\prime}_{2}(1-12J_{2})+\frac{34}{3}J_{3}+\left(\frac{20371}{1728\pi^{2}}-\frac{81}{32\pi^{2}}\ln\frac{m^{2}}{M^{2}}\right)J_{2}\right],$
(47)

therefore still to be solved numerically as an implicit function since
$J_{i}\equiv J_{i}(\frac{m}{T},\frac{\mu}{T})$.

Figure 5: Comparison of NLO RGOPT quark pressure Eq.(26) with
$\overline{m}(g^{2})$ (green, thin lines), LO RGOPT (dotdashed), NLO PT (blue,
dashed), NLO HTLpt quark pressure (red, dotted) with scale dependence $\pi
T\leq M\leq 4\pi T$ (bands) and central scale $M=2\pi T$ (lines) at
$\mu_{B}=0$.

The above expression readily gives a real solution, and allows to consider
$\mu\neq 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.(47) into the RGOPT NLO quark pressure Eq.(26)
(keeping also exact thermal integrals consistently in the latter), gives the
results illustrated for $\mu=0$ in Fig. 5, compared with the standard NLO PT
pressure Eq.(8), and also with the NLO HTLpt (quark) pressure. [NB for a
consistent comparison with the latter at this stage, we have extracted only
the quark contributions within the complete QCD NLO HTLpt pressure, which is
not a trivial separation as in the case of NLO pQCD. Details are explained in
Appendix C].

Alternatively, proceeding similarly with the RG Eq.(43) and (46) gives

$\frac{m^{2}_{RG}}{T^{2}}=\frac{9}{2}gJ_{2}+\frac{g^{2}}{32\pi^{2}}\,\left(172-81\ln\frac{m^{2}}{M^{2}}\right)J_{2},$
(48)

observing that the LO term and the $\ln M$ dependence are identical to those
in Eq.(47). This illustrates that although the MOP and RG prescriptions are
quite different if considering their exact determinations, perturbatively they
differ only by ${\cal O}(g^{2})$ terms, thus formally higher order than the
original NLO perturbative pressure from which they were both constructed.
Moreover, inserting Eq.(48) within Eq.(26) gives almost identical results as
in Fig. 5. Note also that in both Eqs.(47) and (48) the running $g(M)$ exactly
cancels the $M$-dependence at ${\cal O}(g^{2})$, as easily checked using
Eqs.(25), (15), and (52).

As seen in Fig. 5 the RGOPT pressure with the (MOP or RG)
$\overline{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.(47) or Eq.(48).
Actually, it is rather closer to the higher orders standard pQCD pressure, as
will be illustrated below, partly due to Eq.(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 999We mention that the NLO HTLpt pressure in
Fig. 5 (and similarly below in Figs. 6-7, Figs. 9-10) is somewhat different
than the results in Ref.[23], specially at very low $T$. This is due to
considering here only its pure quark contributions, and partly also from using
the exact Eq.(27) instead of a more approximate two-loop running expression
used in [23].. 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 $\mu\neq 0$, that we omit to illustrate. We will see below that
the more elaborate untruncated RGOPT pressure, accounting for higher orders in
$\overline{m}(g)$, has a more drastically improved scale dependence, which is
a main expected RGOPT feature.

###  V.2 Hot quark matter: $T\neq 0$, $\mu=0$

####  V.2.1 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 $\mu=0$) the exact NLO RGOPT pressure
$P(\overline{m},g,T,\mu)$ obtained from our first $\overline{m}_{MOP}$
prescription, defined by solving Eqs.(35),(42) (as explained in details
Subsec.IV.3). In Fig. 6 the pressure is displayed as function of the
temperature, compared with the LO RGOPT and the standard NLO pQCD Eq.(8), for
the scale dependence $\pi T\leq M\leq 4\pi T$. The reduction of scale
dependence stemming from the now exact (untruncated) NLO RGOPT appears
substantial (about a factor $\sim 2$ improvement for e.g. $T\sim 1$ GeV). The
HTLpt NLO (quark) pressure[23] 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[23]), while a more
drastic departure from the ideal gas is only obtained at NNLO for HTLpt[25].

Figure 6:  RGOPT quark pressures as function of temperature at LO and NLO (MOP
prescription) compared with standard NLO PT (pQCD) and NLO HTLpt pressures,
with scale dependence $\pi T\leq M\leq 4\pi T$ at $\mu_{B}=0$.

Figure 7: Same captions as for Fig.6 but with the RGOPT pressure obtained from
alternative $\overline{m}_{RG}$ prescription Eqs.(43), (46).

####  V.2.2 Alternative RG prescription

Similarly to Fig. 6, we illustrate in Fig. 7 the exact NLO RGOPT pressure as
obtained from the alternative $\overline{m}_{RG}$ prescription defined from
solving Eqs.(43) and (46) (explained in details in Subsec.IV.4). As is seen
the RGOPT reduction of remnant scale dependence is even more substantial than
for the previous $\overline{m}_{MOP}$ 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, $\Delta
P/P\equiv(P(M=4\pi T)/P(M=\pi T)-1$, for the different approximations as
indicated.

Figure 8:  $\Delta P/P\equiv P(M=4\pi T)/P(M=\pi T)-1$ as function of
temperature (for $\mu_{B}=0$) for the different NLO RGOPT prescriptions
compared to standard NLO pQCD, with scale dependence $\pi T\leq M\leq 4\pi T$.

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
reasonable crosscheck of the moderate dependence upon the details of the
optimization prescriptions, already observed here at NLO. For both the MOP and
RG prescriptions lower pressure values are obtained at moderate temperatures
as compared to LO RGOPT, NLO HTLpt and NLO pQCD in Figs 6, 7.

###  V.3 Hot and dense quark matter

Figure 9:  RGOPT pressure as function of the temperature at LO and NLO (MOP
prescription), compared with NLO pQCD and NLO HTLpt pressures, with scale
variation $\pi\sqrt{T^{2}+\mu^{2}/\pi^{2}}\leq M\leq
4\pi\sqrt{T^{2}+\mu^{2}/\pi^{2}}$ at $\mu_{B}=1.2$ GeV.

Figure 10: Same captions as in Fig. 9 with alternative NLO RG prescription.

We now consider a nonzero chemical potential values. Since the MOP (35), (42)
and RG (43), (46) prescriptions are defined quite generically they can be
readily applied to the more general $T,\mu\neq 0$ case. As a representative
physical value we illustrate our results for $\mu_{B}=1.2$ GeV. For the
renormalization scale variation range we take as is common
$\pi\sqrt{T^{2}+\mu^{2}/\pi^{2}}\leq M\leq 4\pi\sqrt{T^{2}+\mu^{2}/\pi^{2}}$
within the exact NLO running coupling Eq.(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
$\mu_{B}$ value the qualitative picture is very similar to the $\mu_{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.

###  V.4 Including glue contribution: confrontation to lattice results

In principle a rather similar RGOPT treatment of the pure glue sector should
be possible, generalizing within our framework the hard thermal loop (HTL)
originally proposed in [11], that essentially introduces a gauge-invariant
(non-local) effective Lagrangian properly describing a gluon (thermal) “mass”
term in the HTLpt approximation. However, this requires technically the
evaluation of presently unknown and quite involved thermal integrals, that we
leave for future work [35]. 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 [40]

$\frac{P^{PT}_{g}}{P_{g,SB}}=1-\frac{15}{4}\left(\frac{g}{4\pi^{2}}\right)+{\cal
O}(g^{2}),$ (49)

where the ideal gluon gas pressure is $P_{g,SB}=(8\pi^{2}/45)\,T^{4}$. Thus we
simply add the perturbative NLO contribution Eq.(49) (properly normalized) to
our NLO RGOPT quark contributions Eq.(26), and for the numerical illustrations
below we normalize our results to the full ideal pressure of quarks plus
gluons: $P_{SB}\to P_{q,SB}+P_{g,SB}$ 101010As a slight abuse of notation,
note that in Figs. 5-10 where only quark contributions are included, $P_{SB}$
designates the sole quark ideal pressure Eq.(9), while in Figs.11-13 below
$P_{SB}\equiv P_{q,SB}+P_{g,SB}$.
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 $\overline{m}$,
Eq.(47), for the quark contribution, but supplemented now by the NLO glue
contribution, Eq.(49). The resulting RGOPT pressure is compared with both the
(massless quark) state-of-the-art ${\rm N}^{3}$LO pQCD, which expression is
taken from Ref. [9], and to available LQCD results from Ref. [41, 42, 43]. 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.(48)
gives almost undistinguishable results from Fig. 11, illustrating the low
order perturbative consistency of the two different MOP and RG prescriptions).

Figure 11: RGOPT $P(\overline{m}(g^{2}))$ plus NLO $P^{PT}_{g}$ pressure as
function of $T$ (green band) compared to (${\rm N}^{3}$LO, $g^{3}\ln g$) pQCD
(light blue band), with scale dependence $\pi T\leq M\leq 4\pi T$, and to
lattice data [41, 42, 43] at $\mu_{B}=0$.

Figure 12:  Full NLO RGOPT (MOP prescription) plus NLO $P^{PT}_{g}$ pressure
as function of $T$ (grey band) compared to (${\rm N}^{3}LO\,g^{3}\ln g$) pQCD
(light blue band), with scale dependence $\pi T\leq M\leq 4\pi T$, and to
lattice data [41, 42, 43] at $\mu_{B}=0$.

Figure 13:  Full NLO RGOPT (RG prescription) plus NLO $P^{PT}_{g}$ pressure
(brown band) compared to ${\rm N}^{3}{\rm LO}\,g^{3}\ln g$ pQCD (light blue
band), NLO HTLpt (light green band) and NNLO HTLpt (light red band), with
scale dependence $\pi T\leq M\leq 4\pi T$, and to lattice data [41, 42, 43] at
$\mu_{B}=0$.

Next in Figs.12 and 13, we illustrate similarly the results obtained upon
adding the NLO PT glue contributions Eq.(49) to the NLO RGOPT quark pressure
respectively for the (exact) MOP and RG prescriptions. These are compared with
the state-of-the-art ${\rm N}^{3}$LO pQCD [9], and to LQCD results [41, 42,
43]. 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 [23] and the state-of-the-art NNLO HTLpt [25]. The pressure from the
RG prescription gives the smallest residual scale uncertainties, and is in
remarkable agreement with LQCD data in [41] for the central scale $M=2\pi T$,
for temperatures as low as $T\sim 0.25$ GeV up to $T=1$ GeV, the highest value
considered in [41]. (More precisely let us mention that for the five available
LQCD points in [41] with $T>0.3\,{\rm 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 [42] at intermediate $T$. The RGOPT pressure is somewhat closer to LQCD
results from [41] than the NNLO HTLpt pressure for $0.5\,{\rm GeV}\lesssim
T\lesssim 1\,{\rm GeV}$, while at higher $T$ values HTLpt is nearer to the
results of [43], and RGOPT shows more sizeable differences of order $5-7\%$. A
concomitant feature however is the visible tension between low [41] and higher
$T$ [43] LQCD data in their common temperature range 111111We show LQCD data
as given in publicly available files[41, 42, 43], that do not include
systematic uncertainties..

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.2 and
used in Ref.[34], 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.[34]). In the present $T\neq 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 $\sim g\,m^{4}\,s_{2}$ subtraction term of
Eq.(16), that is formally of NLO ${\cal O}(g)$ 121212 This variant is the next
order analogue of including the NLO coefficient $s_{1}\neq 0$ within LO RGOPT,
see e.g. Eq.(30).. The $s_{2}$ expression [31, 32] incorporates three-loop
order RG coefficient dependence, thus for consistency we took a three-loop
perturbative running coupling generalizing Eq.(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

$\Delta\equiv\varepsilon-3P=T\frac{\partial P}{\partial
T}-4P=T^{5}\,\partial(P/T^{4})/\partial T,$ (50)

(where the second and third equalities are of course valid only for $\mu=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[41, 42, 43] only. A very good
agreement with LQCD results of [41, 42] is obtained for $0.3\,{\rm
GeV}\lesssim T\lesssim 1\,{\rm GeV}$, while there are more visible differences
with the higher $T$ results from [43]. 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.

Figure 14:  NLO RGOPT (RG prescription) trace anomaly
$\Delta\equiv\varepsilon-3P$ (including $\Delta^{PT}_{g}$) (brown band)
compared to lattice data [41, 42, 43]. The additional dashed lines illustrate
the scale uncertainty originating solely from RGOPT quark contributions within
the full scale uncertainty added by $\Delta_{g}^{PT}$ (brown) band.

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[11] (in
particular with a thermal gluon mass parameter, $m_{D}$). There are
essentially three important differences:

  * •

First, the perturbative RG-restoring subtraction terms, like in Eq.(16)
typically, are missing in HTLpt. Accordingly the latter lacks perturbative RG-
invariance formally by a leading order term of the massive theory pressure,
${\cal O}(m^{4})\ln(M/m)$. Now since for any (gluon or quark) thermal masses,
$m^{2}\sim\\#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.[22, 24, 25].

  * •

Second, the interpolating Lagrangian used in HTLpt is linear, namely with an
exponent $a=1$ in the HTL equivalent of Eq.(20), instead of our RG-determined
Eq.(23). As we have shown[26] 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.(21) in
HTLpt (as was done e.g. in [23, 24]), nonreal $\overline{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 $\overline{m}_{D}$ arbitrary mass by a perturbative thermal mass
[22, 25], and taking the quark mass $\overline{m}_{q}=0$. However, enforcing
perturbative masses is partly lacking the a priori more nonperturbative
behaviour rooted in 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
nonperturbative approximations with consistent RG properties already at LO
(one-loop). Our NLO results have been compared to NLO and state-of-the-art
${\rm N}^{3}\mbox{LO}$ pQCD predictions as well as to the state-of-the-art
(NNLO) HTLpt results. Scale variations in the range $\pi T\leq M\leq 4\pi 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, ${\cal 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.(47), (48)), and is substantially reduced within
the resulting physical pressures. Using the RG Eq.(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 [26]). Note that 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 some present technical limitations in
applying the RGOPT plainly to the glue sector, in this work we have adopted a
simple-minded approach, adding the formally same order 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\pi T$ results are
in remarkable agreement with the LQCD results [41, 42] for temperatures as low
as $T\sim 0.25$ GeV, which lies within the nonperturbative regime, up to $T=1$
GeV. The striking matching with LQCD results from Ref. [41] 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\,{\rm
GeV}\lesssim T\lesssim 2\,{\rm GeV}$ LQCD data[43]. We remark that the LQCD
pressure results in [41] and in Ref. [43] appear to be in tension in their
common temperature range, while the trace anomaly shows more continuity 131313
Note that the LQCD simulations in Refs. [41, 42, 43] primarily calculate the
trace anomaly $\Delta$, the pressure being derived by the integral method,
i.e. essentially from numerically integrating the last equality in Eq.(50)., 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
degradation to sizable chemical potential values, $\mu_{B}=1.2$ GeV, 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.

###### Acknowledgements.

We thank Peter Petreczky for bringing the results of Ref. [43] to our
attention. We thank Eduardo Fraga and Rudnei Ramos for related discussions.
M.B.P. is partially supported by Conselho Nacional de Desenvolvimento
Científico e Tecnológico (CNPq-Brazil), Process No. 303846/2017-8, and by
Coordenação de Aperfeiçoamento de Pessoal de Nível Superior-(CAPES-
Brazil)-Finance Code 001. This author also thanks the Charles Coulomb
Laboratory, in Montpellier, for the hospitality. T.E.R. thanks the support and
hospitality of CFisUC where part of this work was developed and acknowledges
Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq-Brazil)
and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES-Brazil)
for PhD grants at different periods of time. This work was financed in part by
INCT-FNA (Process No. 464898/2014-5).

##  Appendix A high-$T$ limit

We give here for completeness the well-known $T\gg m,\mu=0$ approximations
(see e.g.[7, 36]) of the basic thermal integrals defined in Eqs (2)-(4):

$2J_{1}(T\gg
m,\mu=0)\approx\frac{7\pi^{2}}{180}-\frac{m^{2}}{12\,T^{2}}+\frac{m^{4}}{T^{4}}\frac{2}{(4\pi)^{2}}\left[\frac{3}{4}-\ln\left(\frac{me^{\gamma_{E}}}{\pi
T}\right)\right]+{\cal O}\left(\frac{m^{6}}{T^{6}}\right)\;,$ (51) $J_{2}(T\gg
m,\mu=0)\approx\frac{1}{12}+\frac{1}{4\pi^{2}}\frac{m^{2}}{T^{2}}\left[\ln\left(\frac{me^{\gamma_{E}}}{\pi
T}\right)-\frac{1}{2}\right]+{\cal O}\left(\frac{m^{4}}{T^{4}}\right)\;,$ (52)

and the more complicated genuine two-loop integral $J_{3}$ of Eq.(4) has a
finite $m\to 0$ limit (however not analytically integrable to our knowledge,
we give below its numerically integrated approximate value):

$J_{3}\left(\frac{m}{T}\to
0,\frac{\mu}{T}=0\right)=\frac{4}{(2\pi)^{4}}\int_{0}^{\infty}d\hat{p}\int_{0}^{\infty}d\hat{q}\,n_{F}(\hat{p})n_{F}(\hat{q})\,\ln\left(\frac{|\hat{p}-\hat{q}|}{\hat{p}+\hat{q}}\right)+{\cal
O}\left(\frac{m^{2}}{T^{2}}\right)\simeq-0.00129532+{\cal
O}\left(\frac{m^{2}}{T^{2}}\right)$ (53)

where $\hat{p},\hat{q}\equiv p/T,q/T$ and $n_{F}(p)=(e^{p}+1)^{-1}$ is the
Fermi-Dirac distribution.

##  Appendix B Numerical $\overline{m}$ solutions at NLO

We discuss here in some details the behavior of the exact NLO numerical
solutions for the two MOP or RG prescriptions as defined in Secs.IV.3, and
IV.4. Note that using directly the MOP Eq.(35) or the RG Eq.(43) makes the AF
solution identification obvious. Concerning the MOP Eq.(35), once $B_{2}$ is
consistently determined by Eq.(42) such as to recover $D_{mop}>0$ in Eq.(37),
one sees from the structure of (35) that $(-)$ (AF) solutions only exist if
$-\ln(m^{2}/M^{2})+B_{mop}>0$, and conversely $(+)$ (non-AF) solutions only
exist if $-\ln(m^{2}/M^{2})+B_{mop}<0$.

Figure 15:  AF and non-AF roots of MOP Eqs.(35), (42) for $M=2\pi T$ and for
two representative $T$ values, $T=0.5$ GeV (dashed), $T=1$ GeV (thick).

Once $M,g(M)$ are taken to be $T,\mu$-dependent via the perturbative running
coupling Eq.(27), Eq.(35) becomes a function of $m/T$ and
$g(T/\Lambda_{\overline{MS}})$. Despite the nonlinear dependence in $m/T$, at
the level of Eq.(35) both the AF and non-AF solutions happen to be unique in
their respective existence range. This is illustrated in Fig. 15 (for $\mu=0$)
for two representative low to moderate temperatures, respectively $T=0.5$ and
$T=1$ GeV, and for the central scale choice $M=2\pi T$. It is also clear that
for any $T$ the smallest solution is the AF one: Indeed for $g(\pi T\leq M\leq
4\pi 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
$\overline{m}(\mbox{AF})<m_{0}<\overline{m}(\mbox{non-AF})$. The value of
$m_{0}$ depends quite strongly on $T$ (and $M$): typically for the input
corresponding to Fig. 15 with $M=2\pi T$, one finds $m_{0}\simeq 1.28(1.91)$
for $T=0.5(1.0)$ GeV respectively. (Note also that in Fig. 15 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 $\mu=0$, following the AF-matching $\overline{m}$ of Eq.(35) 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\sim\pi T$ one hits on $M\sim\Lambda_{\overline{MS}}$ where the
perturbative coupling diverges. For sizable $\mu\neq 0$ the latter problem if
avoided if defining as conventional $M\sim\pi\sqrt{T^{2}+\mu^{2}/\pi^{2}}$
(provided that one is not in the case of both $T\ll\mu$ and small $\mu$).

Finally concerning the RG Eq.(43), both NLO solutions are already AF-matching,
giving thus a unique solution upon using the prescription Eq.(46). Numerically
the exact $\overline{m}_{RG}$ solution of Eq.(43) is somewhat larger than
$\overline{m}_{MOP}$ for a given $T$, as illustrated in Fig. 3.

##  Appendix C NLO and NNLO HTLpt expressions

For completeness we specify here how the NLO[23] and NNLO[25] 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 [23] to the pure glue NLO HTLpt pressure (given e.g. in
Eq.(4.8) of second Ref. in [22]), 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 [23], that amounts
to remove in these expressions the pure glue contributions (terms $\propto
c_{A}$ in Eq.(55) in [23]).
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 Ref.[25] (see also Ref.[24]).

## References

  * [1] P. de Forcrand, PoS LAT 2009, 010 (2009); G. Aarts, J. Phys. Conf. Ser. 706, 022004 (2016).
  * [2] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz and K. K. Szabo, Nature 443, 675 (2006); Y. Aoki, S. Borsanyi, S. Durr, Z. Fodor, S. D. Katz, S. Krieg and K. K. Szabo, JHEP 06, 088 (2009); S. Borsanyi et al. [Wuppertal-Budapest], JHEP 09, 073 (2010); A. Bazavov, T. Bhattacharya, M. Cheng, C. DeTar, H. T. Ding, S. Gottlieb, R. Gupta, P. Hegde, U. M. Heller and F. Karsch, et al. Phys. Rev. D 85, 054503 (2012).
  * [3] R. Machleidt and D. R. Entem, Phys. Rept. 503, 1 (2011).
  * [4] Y. Nambu and G. Jona-Lasinio, Phys. Rev. 122, 345 (1961).
  * [5] M. Buballa, Phys. Rept. 407, 205 (2005).
  * [6] J. P. Blaizot, E. Iancu and A. Rebhan, In *Hwa, R.C. (ed.) et al.: Quark gluon plasma* 60-122 [hep-ph/0303185]; U. Kraemmer and A. Rebhan, Rept. Prog. Phys. 67, 351 (2004).
  * [7] M. Laine and A. Vuorinen, Lect. Notes Phys. 925, 1 (2016).
  * [8] J. Ghiglieri, A. Kurkela, M. Strickland and A. Vuorinen, Phys. Rept. 880, 1 (2020).
  * [9] K. Kajantie, M. Laine, K. Rummukainen and Y. Schroder, Phys. Rev. D 67, 105008 (2003).
  * [10] A. Vuorinen, Phys. Rev. D 68, 054017 (2003).
  * [11] E. Braaten and R. D. Pisarski, Phys. Rev. D 45, 1827 (1992).
  * [12] F. Karsch, A. Patkos and P. Petreczky, Phys. Lett. B 401, 69 (1997); S. Chiku and T. Hatsuda, Phys. Rev. D 58, 076001 (1998); J. O. Andersen, E. Braaten and M. Strickland, Phys. Rev. D 63, 105008 (2001); J. O. Andersen and M. Strickland, Phys. Rev. D 64, 105012 (2001); J. O. Andersen and M. Strickland, Annals Phys. 317, 281 (2005).
  * [13] J. O. Andersen, E. Braaten and M. Strickland, Phys. Rev. Lett. 83, 2139 (1999); J. O. Andersen, E. Braaten and M. Strickland, Phys. Rev. D 61, 074016 (2000).
  * [14] P.M. Stevenson, Phys. Rev. D 23, 2916 (1981); Nucl. Phys. B 203, 472 (1982).
  * [15] A. Okopinska, Phys. Rev. D 35, 1835-1847 (1987) doi:10.1103/PhysRevD.35.1835;
  * [16] H. Yamada, Z. Phys. C 59, 67-76 (1993) doi:10.1007/BF01555840
  * [17] A. Duncan and M. Moshe, Phys. Lett. B 215, 352-358 (1988) doi:10.1016/0370-2693(88)91447-5; H. F. Jones and M. Moshe, Phys. Lett. B 234, 492-496 (1990) doi:10.1016/0370-2693(90)92045-K .
  * [18] R. P. Feynman and H. Kleinert, Phys. Rev. A 34, 5080 (1986); H. Kleinert, Phys. Rev. D 57, 2264 (1998); Phys. Lett. B 434, 74 (1998); Phys. Rev. D 60, 085001 (1999); Mod. Phys. Lett. B 17, 1011 (2003).
  * [19] J. O. Andersen, L. Kyllingstad and L. E. Leganger, JHEP 0908, 066 (2009).
  * [20] S. K. Gandhi, H. F. Jones and M. B. Pinto, Nucl. Phys. B 359, 429 (1991).
  * [21] J. L. Kneur, M. B. Pinto and R. O. Ramos, Phys. Rev. D 74, 125020 (2006).
  * [22] J. O. Andersen, M. Strickland and N. Su, Phys. Rev. Lett. 104, 122003 (2010); JHEP 1008, 113 (2010).
  * [23] N. Haque, M. G. Mustafa and M. Strickland, Phys. Rev. D 87, 105007 (2013).
  * [24] J. O. Andersen, L. E. Leganger, M. Strickland and N. Su, JHEP 1108, 053 (2011); S. Mogliacci, J. O. Andersen, M. Strickland, N. Su and A. Vuorinen, JHEP 1312, 055 (2013); N. Haque, J. O. Andersen, M. G. Mustafa, M. Strickland and N. Su, Phys. Rev. D 89, 061701 (2014).
  * [25] N. Haque, A. Bandyopadhyay, J. O. Andersen, M. G. Mustafa, M. Strickland and N. Su, JHEP 1405, 027 (2014).
  * [26] J. L. Kneur and M. B. Pinto, Phys. Rev. D 92, 116008 (2015).
  * [27] J. L. Kneur and M. B. Pinto, Phys. Rev. Lett. 116, 031601 (2016).
  * [28] J. L. Kneur and A. Neveu, Phys. Rev. D 81, 125012 (2010).
  * [29] J. L. Kneur and A. Neveu, Phys. Rev. D 88, 074025 (2013).
  * [30] M. Tanabashi et al. [Particle Data Group], Phys. Rev. D 98, 030001 (2018).
  * [31] J. L. Kneur and A. Neveu, Phys. Rev. D 92, 074027 (2015).
  * [32] J. L. Kneur and A. Neveu, Phys. Rev. D 101, 074009 (2020).
  * [33] G. N. Ferrari, J. L. Kneur, M. B. Pinto and R. O. Ramos, Phys. Rev. D 96, 116009 (2017).
  * [34] J. L. Kneur, M. B. Pinto and T. E. Restrepo, Phys. Rev. D 100, 114006 (2019).
  * [35] J. L. Kneur and M. B. Pinto; in preparation.
  * [36] J. I. Kapusta and C. Gale, “Finite-temperature field theory: Principles and applications” (Cambridge University Press, 2006).
  * [37] M. Laine and Y. Schröder, Phys. Rev. D 73, 085009 (2006).
  * [38] J. L. Kneur, Phys. Rev. D 57, 2785 (1998).
  * [39] I. T. Drummond, R. R. Horgan, P. V. Landshoff and A. Rebhan, Nucl. Phys. B 524, 579 (1998).
  * [40] E. V. Shuryak, Sov. Phys. JETP 47, 212 (1978); S. A. Chin, Phys. Lett. B 78, 552 (1978).
  * [41] S. Borsanyi, G. Endrodi, Z. Fodor, A. Jakovac, S. D. Katz, S. Krieg, C. Ratti and K. K. Szabo, JHEP 11, 077 (2010).
  * [42] S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg and K. K. Szabo, Phys. Lett. B 730, 99 (2014).
  * [43] A. Bazavov, P. Petreczky and J. H. Weber, Phys. Rev. D 97, 014510 (2018).