# Renormalization group improved pressure for hot and dense quark matter

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

► To cite this version:

Jean-Loïc Kneur, Marcus Benghi Pinto, Tulio E. Restrepo. Renormalization group improved pressure for hot and dense quark matter. Physical Review D, 2021, 104 (3), pp.034003. 10.1103/Phys-RevD.104.034003. hal-03129450

HAL Id: hal-03129450

https://hal.science/hal-03129450

Submitted on 1 Jun 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Renormalization group improved pressure for hot and dense quark matter

Jean-Loïc Kneur<sup>*</sup>,<sup>1,∗</sup> Marcus Benghi Pinto<sup>2,†</sup> and Tulio E. Restrepo<sup>2,‡</sup>

<sup>1</sup>Laboratoire Charles Coulomb (L2C), UMR 5221 CNRS-Université Montpellier,
34095 Montpellier, France

<sup>2</sup>Departamento de Física, Universidade Federal de Santa Catarina, 88040-900 Florianópolis, SC, Brazil

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

(Received 21 January 2021; accepted 6 July 2021; published 9 August 2021)

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 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 \lesssim T \lesssim 1$ GeV, with a remnant scale dependence drastically reduced
as compared to HTLpt.

DOI: 10.1103/PhysRevD.104.034003

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 Tpc ≈ 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 nonper-
turbative 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 by-product,
one should observe a critical end point (CP) signaled 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 approxi-
mation (which is related to MFT) may produce misleading
results at criticality. More analytical alternatives to LQCD
can partly address the deconfinement and/or chiral sym-
metry 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 incorpo-
rating more basic QCD dynamics in well-defined approx-
imations, 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 temper-
ature 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 temper-
ature, 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

Published by the American Physical Society under the terms of
the Creative Commons Attribution 4.0 International license.
Further distribution of this work must maintain attribution to
the author(s) and the published article's title, journal citation,
and DOI. Funded by SCOAP³.

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

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

<footnote>tulio.restrepo@posgrad.ufsc.br</footnote>
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–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 occurring in Yang-Mills theory, intrinsically limiting in practice the weak coupling expansion to maximal order $\alpha_s^3 \ln \alpha_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–24] (as we dub it here), linear $\delta$ 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–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 \gtrsim 2T_{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 [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–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 $\alpha_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 $\phi^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. Section VI contains our conclusions and perspectives. Finally, three appendixes 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$,<fnref n="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 = \mu = 0$) results of Ref. [39] and the $T, \mu \neq 0$ results of Refs. [44,45]. Considering the case of degenerate masses $m_u = m_d = m_s \equiv m$, the renormalized pressure in the $\overline{\text{MS}}$ renormalization scheme, normalized per flavor, is

$$ \frac{P_1^{PT}}{N_f N_c} = -\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(2\pi)^4} C_F \left(L_m^2 - \frac{4}{3}L_m + \frac{3}{4}\right) \\ -g C_F \left\{ \left[ \frac{m^2}{4\pi^2} (2-3L_m) + \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\}, \quad (2.1) $$

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 \hat{\mathbf{p}}}{(2\pi)^3} \left\{ \ln[1 + e^{-(E_p + \frac{m}{T})}] + \ln[1 + e^{-(E_p - \frac{m}{T})}] \right\}, \quad (2.2) $$

$$ \text{with } \hat{\mathbf{p}} \equiv \mathbf{p}/T, \quad E_p = \sqrt{\hat{\mathbf{p}}^2 + \frac{m^2}{T^2}}, $$

$$ P^{PT} = \text{[Diagram]} + \text{[Diagram with self-coupling]} + \mathcal{O}(g^2) $$

<caption>FIG. 1. Feynman diagrams contributing to the perturbative quark pressure up to NLO $\mathcal{O}(g)$.</caption>

$$ J_2\left(\frac{m}{T}, \frac{\mu}{T}\right) = \int \frac{d^3 \hat{\mathbf{p}}}{(2\pi)^3} \frac{1}{E_p} [f^+(E_p) + f^-(E_p)], \quad (2.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_0^\infty \int_0^\infty \frac{d\hat{p} \, \hat{p} \, d\hat{q} \, \hat{q}}{E_p E_q} \\ \times \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] \right. \\ \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] \right\}, \quad (2.4) $$

where

$$ \Sigma_{\pm} = f^{+}(E_{p})f^{\pm}(E_{q}) + f^{-}(E_{p})f^{\mp}(E_{q}), \quad (2.5) $$

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

$$ f^{\pm}(E_p) = \frac{1}{1 + e^{(E_p \pm \frac{m}{T})}}, \quad (2.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_1^{PT}(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) \quad (2.8) $$

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

<footnote n="1">¹In all what follows we normalize for convenience the (running) coupling in the $\overline{\text{MS}}$-scheme as $g(M) \equiv 4\pi\alpha_S(M)$.</footnote>
$$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). \quad (2.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, \quad (2.10)$$

where our normalization convention for the QCD $\beta$-function and anomalous mass dimension $\gamma_m$ is

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

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

where to two-loop order,

$$(4\pi)^2 b_0 = 11 - \frac{2}{3} N_f, \quad (2.13)$$

$$(4\pi)^4 b_1 = 102 - \frac{38}{3} N_f, \quad (2.14)$$

$$\gamma_0 = \frac{1}{2\pi^2}, \quad (4\pi)^4 \gamma_1 = \frac{404}{3} - \frac{40}{9} N_f. \quad (2.15)$$

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,

$$\frac{P^{\text{RGPT}}}{N_c N_f} = \frac{P^{\text{PT}}}{N_c N_f} - \frac{m^4}{g} \sum_k s_k g^k, \quad (2.16)$$

where the $s_i$ are determined at successive perturbative orders so that

$$M \frac{d}{dM} \left( \frac{P^{\text{RGPT}}}{N_c N_f} \right) = \mathcal{O}(g^2 m^4), \quad (2.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. (2.10). One finds explicitly [39,42]

$$s_0 = -[(4\pi)^2(b_0 - 2\gamma_0)]^{-1}, \quad (2.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]. \quad (2.19)$$

## 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^{\text{RGPT}}(m, g)$ in Eqs. (2.16), (2.18), (2.19). The subtraction contributions can be viewed as extending the Lagrangian by vacuum terms $\sim(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 $\sim 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

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

$m$ being now an arbitrary mass. Equation (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^{\text{RGPT}}(m, g)$ of Eq. (2.16), thus incorporating the previously obtained vacuum subtractions.<fnref n="2" />

(3) The resulting expression from Eq. (2.20) is expanded in powers of $\delta$, at the same perturbative order considered, and next setting $\delta \to 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),

$$\left. \frac{\partial}{\partial m} P^{\text{RGOPT}}(\mathcal{O}(\delta^k), \delta \to 1) \right|_{\bar{m}} = 0. \quad (2.21)$$

Equation (2.21) is however not the sole nor the most compelling prescription once RG invariance properties are considered. Moreover, while in related OPT/ SPT or HTLpt approaches the simple (linear) value $a = 1$ was used in Eq. (2.20) mainly for simplicity, it was observed [35,37] that step 3 generally spoils RG invariance at leading order, even when Eq. (2.20) is

<footnote n="2">²For overall RG consistency it is crucial to perform first the perturbative RG-restoring procedure as in Eq. (2.16), before modifying the perturbative expansion according to Eq. (2.20), as the two procedures do not commute.</footnote>
performed on a perturbatively RG invariant pressure like Eq. (2.16).

(4) Thus in conjunction with step 3, to preserve good RG properties of the variationally modified pressure the further RGOPT ingredient is to fix [35,37] the exponent a purposely introduced in Eq. (2.20): we require the $\delta \to 1$ pressure to satisfy, since Eq. (2.21) is used, the reduced (massless) RG equation,

$$[M\partial_M + \beta(g)\partial_g]P_{\delta \to 1}^{\text{RGOPT}} = 0. \quad (2.22)$$

For this it is sufficient to expand to leading order ($\delta^0$), yielding

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

that only depends on the universal (scheme-independent) LO RG coefficients,<fnref n="3" /> in agreement with previous RGOPT applications to which we refer for detailed demonstration [34,37,42]. At NLO Eq. (2.22) is no longer exactly satisfied by Eq. (2.23) alone, thus it can provide an alternative determination of $\bar{m}$ besides Eq. (2.21), as will be illustrated in concrete NLO applications below. Importantly, keeping Eq. (2.23) at higher orders further guarantees that the only acceptable solutions are those matching [37] the perturbative asymptotic freedom (AF) behavior for $g \to 0$ at $T = 0$. This simple but compelling criterion often selects a unique solution, even at five-loop 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 $\delta^0$ the resulting RGOPT pressure is given, keeping all terms of formally one-loop order, by

$$\begin{aligned} \frac{P_0^{\text{RGOPT}}}{N_f N_c} = & - \frac{2m^4}{(4\pi)^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}{(4\pi)^2 g b_0}. \end{aligned} \quad (2.24)$$

Remark that the LO $s_0$ coefficient, Eq. (2.18), has produced the last term $\propto 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 \neq 0$ as explained. Yet, since $s_1$ enters the pressure formally at $\mathcal{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 \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(M) = \frac{1}{2b_0 \ln(M/\Lambda_{\overline{\text{MS}}})}. \quad (2.25)$$

Proceeding similarly at the next RGOPT order, the NLO pressure reads, after performing steps 1–4 above,

$$\begin{aligned} \frac{P_1^{\text{RGOPT}}}{N_f N_c} = & -\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}{(2\pi)^2} \left(\frac{\gamma_0}{b_0}\right) \left(\frac{1}{2} - L_m\right) + m^2 \left(\frac{\gamma_0}{b_0}\right) T^2 J_2\left(\frac{m}{T}, \frac{\mu}{T}\right) \\ & + \frac{m^4}{(4\pi)^2 b_0} \left\{ \frac{1}{g} \left(1 - \frac{\gamma_0}{b_0}\right) + \left[ (b_1 - 2\gamma_1)\pi^2 - \frac{b_0 - 2\gamma_0}{3} \right] \right\} - 3g C_F \frac{m^4}{2(2\pi)^4} \left( L_m^2 - \frac{4}{3}L_m + \frac{3}{4} \right) \\ & - g C_F \left\{ \frac{m^2}{4\pi^2} (2 - 3L_m) + \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) \end{aligned} \quad (2.26)$$

<footnote n="3">³At NLO and higher orders one could possibly further generalize the interpolation $(1 - \delta)^a$ in Eq. (2.20) with $\delta^2$ and higher order terms without spoiling the crucial LO RG properties entirely fixed by Eq. (2.23). But this would involve extra arbitrary variational parameters with no compelling reasons, with no way of fixing those so that the remnant NLO scale dependence of $\mathcal{O}(g^2 m^4)$ would exactly disappear. Keeping a minimal form with only $a$ in Eq. (2.23) gives the same prescription at successive orders, thus more sensible for comparisons.</footnote>
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]),

$$ \ln \frac{M}{\Lambda_{\overline{\text{MS}}}} = \frac{1}{2b_0 g} + \frac{b_1}{2b_0^2} \ln \left( \frac{b_0 g}{1 + \frac{b_1}{b_0} g} \right), \quad (2.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 [38], $\Lambda_{\overline{\text{MS}}} = 335$ MeV for $N_f = 3$, that equivalently corresponds to $\alpha_s(N_f = 3, 1.5 \text{ GeV}) \approx 0.326$. (NB the latter $\alpha_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 ($\delta^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,

$$ \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, \quad (3.1) $$

whose nontrivial solution gives an RG invariant dressed mass $\bar{m}(g, T, \mu)$, 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 \neq 0$ as given by Eq. (2.19)]. Once inserting $\bar{m}$ in Eq. (2.24) it produces a (one-loop) exactly RG invariant pressure, that takes the compact form,

$$ \frac{P_{0}^{\text{RGOPT}}}{N_f N_c} = 2T^4 J_1\left(\frac{\bar{m}}{T}, \frac{\mu}{T}\right) + \frac{T^2}{2}\bar{m}^2 J_2\left(\frac{\bar{m}}{T}, \frac{\mu}{T}\right) - \frac{\bar{m}^4}{32\pi^2}, \quad (3.2) $$

where it is understood that $\bar{m}$ is the nontrivial solution of Eq. (3.1).<fnref n="4" /> Some properties of the dressed mass $\bar{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 [14] $J_1, J_2$, given in Appendix A. This gives from Eq. (3.1),

$$ \begin{aligned} \bar{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} \\ &\approx \frac{3}{8} g T^2 + \mathcal{O}(g^2), \end{aligned} \quad (3.3) $$

or, equivalently using Eq. (2.25),

$$ \bar{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{\text{MS}}}} \right) \right]^{-1}, \quad (3.4) $$

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

$$ m_{\text{PT}}^2 = \frac{g}{6} T^2 + \mathcal{O}(g^2). \quad (3.5) $$

In contrast $\bar{m}$ in Eq. (3.3) represents an intermediate variational quantity, whose meaning is merely once being inserted in $P(\bar{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. (2.16), leading to $\bar{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.<fnref n="5" /> Accordingly at LO and in the high-$T$ approximation, using Eq. (2.25), Eq. (3.2) takes the simpler form,

$$ \frac{P_{0}^{\text{RGOPT}}}{P_{\text{SB}}} \approx 1 - \frac{5}{14} \left[ \ln \left( \frac{\pi T}{e^{\gamma_E} - \frac{53}{84} \Lambda_{\overline{\text{MS}}}} \right) \right]^{-1}, \quad (3.6) $$

normalized to the SB ideal quark gas $P_{\text{SB}}$ Eq. (2.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. (3.6): at high-$T$ the explicit $M$-dependence in Eq. (3.3) has been automatically traded for a dependence in $g(\sim \pi T / \Lambda_{\overline{\text{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 (nonre-summed) thermal perturbative expansions.

<footnote n="4">⁴Notice that the explicit dependence upon $s_1$ cancelled in $P_{0}^{\text{RGOPT}}$ Eq. (3.2) upon using Eq. (3.1), but the solution $\bar{m}$ of Eq. (3.1) does depend on $s_1$ as specified.</footnote>

<footnote n="5">⁵In the simpler $\mathcal{O}(N)\phi^4$ model, the analogous LO RGOPT [34] resums all large-$N$ contributions, reproducing the exactly known large-$N$ pressure [47], including nonanalytic terms $\sim \lambda^{3p/2}$, $p \ge 1$, typical of a boson gas pressure.</footnote>
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 $\bar{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 $\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.

## 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 $\bar{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. (2.27) and scale $M$ chosen as a combination of $\pi T$ and $\mu$ when both $T, \mu$ are nonzero. As shown above in Sec. III A, the LO RGOPT pressure exhibits one-loop-exact scale invariance as a consequence of the simple form of Eq. (2.24) only depending on $b_0$. While at NLO or beyond, $P^{\text{RGOPT}}(\bar{m}(g), M, \dots)$ 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, $\sim 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 $\gamma_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 $\sim 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 \sim 2\pi T$ (for $\mu = 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(\pi 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 \to 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 $\bar{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{\text{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 = \mu = 0$ case [37] as well as within the $T = 0, \mu \neq 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, \bar{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 [42], it becomes numerically quite challenging for the rather involved $T, \mu \neq 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.

# IV. NLO PRESCRIPTIONS

## 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,

$$m \to m'(1 + B_2 g^2), \qquad (4.1)$$

where a single $B_2$ coefficient parametrizes a perturbative NLO scheme change from the original $\overline{\text{MS}}$-scheme.<fnref n="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. (4.1) consistently on the original PT expression (2.1) prior to its modification induced from Eq. (2.20) with the subsequent $\delta$-expansion. The net RSC modification to the pressure is to add an extra term, $-4g(m')^4 s_0 B_2$, entering thus the resulting exact NLO Eq. (2.21) or Eq. (2.22). Thus Eq. (4.1) modifies the latter equations purposefully, now considering those equations as constraints for the arbitrary mass $m'$, after the modifications from Eq. (2.20).<fnref n="7" /> 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.

## B. AF-compatible NLO dressed mass solutions

To identify some relevant properties of the sought dressed $\bar{m}(g, T, \mu)$ solutions we consider first the MOP Eq. (2.21) more explicitly at NLO, thus applied to Eq. (2.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_{\text{mop}} \mp \frac{2\pi}{3g} \sqrt{D_{\text{mop}}} = 0, \qquad (4.2)$$

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

$$B_{\text{mop}} = -\frac{7\pi^2}{9g} + \frac{5}{6} + 4\pi^2 \left( J'_2 + \frac{T^2}{m^2} J_2 \right), \qquad (4.3)$$

$$D_{\text{mop}} = 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) \\ + 36\pi^2 g^2\left(J_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'\right) \\ + 8\pi^2 g^2 \frac{T^2}{m^2}(3J_2 - 1)J_2' - 48\pi^2 g^2\left(J_3' + \frac{T^2}{m^2}J_3\right), \qquad (4.4)$$

where $J_i \equiv J_i(m^2/T^2, \mu/T)$, $J'_i \equiv \partial_x J_i(x)$, (note that here $x \equiv m^2/T^2$). In Eq. (4.4) we explicitly separated the $T, \mu$-independent part within $D_{\text{mop}}$ in the very first line to make its $T, \mu \to 0$ limit clear [remark also that $D_{\text{mop}}(T=0)$ does not depend on $m$].

One first property of Eq. (4.2) is exhibited from expanding it perturbatively to the first few terms. That gives

$$\ln \frac{\bar{m}^2}{M^2} (-) \approx -\frac{16\pi^2}{9g} + \frac{139}{54} + 8\pi^2 \frac{T^2}{\bar{m}^2} J_2 + \mathcal{O}(g), \qquad (4.5)$$

and

$$\ln \frac{\bar{m}^2}{M^2} (+) \approx \frac{2\pi^2}{9g} - \frac{49}{54} + 8\pi^2 J_2' + \mathcal{O}(g). \qquad (4.6)$$

One easily recognizes that, for $T \to 0$ the leading term for $g \to 0$ in $\bar{m}^2(-)$ has the correct AF behavior: $\ln \frac{\bar{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 Eq. (4.2) 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. (A2)]. After straightforward algebra one obtains, for the first few perturbative expansion terms,

$$\frac{\bar{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) + \mathcal{O}(g^3), \qquad (4.7)$$

where we defined for short,

$$L_T \equiv \ln \left( \frac{Me^{\gamma_E}}{\pi T} \right). \qquad (4.8)$$

As seen the AF-compatible solution $\bar{m}(-)$ has a typical perturbative thermal screening mass behavior $m \sim \sqrt{gT}$, 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

<footnote n="6">⁶Equation (4.1) 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. [37] for details).</footnote>

<footnote n="7">⁷To avoid excessive notation proliferation, in what follows once having performed the replacement implied by Eq. (4.1) we simply rename $m' \to m$ the variational mass to be determined from Eq. (2.21) or Eq. (2.22).</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 2. Perturbatively reexpanded NLO RGOPT pressure $P(T, \mu = 0)$ (red band) compared with standard perturbative NLO pressure Eq. (2.8) (blue band), with scale dependence $\pi T \le M \le 4\pi T$.</figcaption></figure>

Eq. (4.2) with (+) has $\bar{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 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 $\mu \neq 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 $\bar{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 $\overline{MS}$ scheme, i.e., $B_2 \equiv 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) $\bar{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. (4.4) is positive, this loss of real $\bar{m}$ solutions arises solely when considering the exact Eq. (4.2): now since all our 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 $\bar{m}(g, T, \mu)$, recovering real solutions for practically most relevant $g$ values. Note that simply requiring $D_{mop} \ge 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$, $\mu$-independent negative term within Eq. (4.4), $-(47/6) g$. Explicitly that gives

$$B_2 = -\frac{329}{1728\pi^2 g} \qquad (4.9)$$

The latter $B_2$ prescription is very simple, and the resulting $\bar{m}_{MOP}$ solution remains real for practically all physically relevant $g(T, \mu)$ values, while still including nontrivial

<footnote>$^8$In particular the exact NLO $\bar{m}$ as obtained below can be such that $\bar{m}/T > 1$ at sufficiently low $T$ (see Fig. 3), somewhat invalidating the high-$T$ approximation.</footnote>

<footnote>$^9$For example for $T = \mu = 0$ where only the first three terms of Eq. (4.4) are nonvanishing, $D_{mop} \le 0$ for rather moderate $g \ge 2.64$, i.e., $\alpha_s \ge 0.21$.</footnote>
<figure><img src="image_4.png" /><figcaption>FIG. 3. Exact LO RGOPT thermal mass (dot-dashed) compared with exact MOP and RG NLO thermal mass for $\pi \le M \le 4\pi T$ at $\mu_B = 0$.</figcaption></figure>

higher order corrections induced from all remnant terms of Eq. (4.2). Other slightly different $B_2$-fixing prescriptions are possible for $T, \mu \ne 0$, but a notable property is that for different $B_2$ choices, that imply different exact $\bar{m}(B_2)$ solutions, the resulting physical pressure $P(\bar{m}(B_2), B_2, \dots)$ Eq. (2.26) happens to be largely insensitive to those unphysical $B_2$ parameter choices provided that $\bar{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. (4.7): $\bar{m}^2 \sim (3/8)gT^2 + O(g^2)$: as easily checked, $B_2$ from Eq. (4.1) only appears at higher order $O(g^3)$ both in the perturbatively expanded $\bar{m}$ Eq. (4.7) and corresponding reexpanded pressure from Eq. (2.26). In other words once $B_2$ is adjusted to recover a real $\bar{m}$ solution of Eq. (4.2), 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. (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_2g^2$ in Eq. (4.1) resulting from Eq. (4.9) as function of $T$. As an important cross-check, it shows that the departure from the original $\overline{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

<figure><img src="image_5.png" /><figcaption>FIG. 4. RSC parameter $B_2g^2(M)$ for the MOP and RG prescriptions for $\pi T \le M \le 4\pi T$ at $\mu_B = 0$.</figcaption></figure>

$\bar{m}(g, T, \mu)$. Once expressed for $\ln(m^2/M^2)$ it takes a similar quadratic form as Eq. (4.2), conveniently normalized as

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

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, \quad (4.11) $$

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. \quad (4.12) $$

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}$,

$$ D_{rg}(B_2) \equiv 0. \quad (4.13) $$

Equation (4.13) determines $B_2$ trivially using Eq. (4.12), leading to a single real AF-compatible solution $\bar{m}_{RG}$ determined from the first two terms of Eq. (4.10), the latter being still an implicit equation in $m$ for $T, \mu \ne 0$ via $J_2$ entering Eq. (4.11). Equation (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 $\bar{m}(B_2)$ is an intermediate variational parameter without much physical meaning outside its use in the pressure. Here the resulting
$\bar{m}_{RG}(B_2)$ still involves arbitrary higher order contributions, as well as nontrivial $T, \mu$ 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, \mu \neq 0$ are not strongly dependent upon those choices.

The dressed exact thermal mass $\bar{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 $\bar{m}_{MOP}$ 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_2g^2$ obtained from Eq. (4.13) is illustrated in Fig. 4 as function of $T$, and compared to the similar MOP $B_2g^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 \lesssim 0.5$ GeV: as already mentioned above the important feature is that the induced departure from the original $\overline{\text{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, \mu$ dependence, is to expand $\bar{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,<fnref n="10" />

<figure><img src="image_6.png" /><figcaption>FIG. 5. Comparison of NLO RGOPT quark pressure Eq. (2.26) with $\bar{m}(g^2)$ (green, thin lines), LO RGOPT (dotdashed), NLO PT (blue, dashed), NLO HTLpt quark pressure (red, dotted) with scale dependence $\pi T \le M \le 4\pi T$ (bands) and central scale $M = 2\pi T$ (lines) at $\mu_B = 0$.</figcaption></figure>

$$ \frac{m_{\text{MOP}}^2}{T^2} = \frac{9}{2} g J_2 + g^2 \left[ \frac{17}{9} J'_2 (1 - 12 J_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], \quad (5.1) $$

therefore still to be solved numerically as an implicit function since $J_i \equiv J_i(\frac{m}{T}, \frac{\mu}{T})$. The above expression readily gives a real solution, and allows us 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. (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 $\mu = 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).

Alternatively, proceeding similarly with the RG Eqs. (4.10) and (4.13) gives

$$ \frac{m_{\text{RG}}^2}{T^2} = \frac{9}{2} g J_2 + \frac{g^2}{32\pi^2} \left( 172 - 81 \ln \frac{m^2}{M^2} \right) J_2, \quad (5.2) $$

observing that the LO term and the $\ln M$ dependence are identical to those in Eq. (5.1). This illustrates that although the MOP and RG prescriptions are quite different if considering their exact determinations, perturbatively they differ only by $O(g^2)$ terms, thus formally higher order than the original NLO perturbative pressure from which they were both constructed. Moreover, inserting Eq. (5.2) within

<footnote n="10">¹⁰As an algebraic subtlety, one should first expand perturbatively the (AF-matching) Eq. (4.2) with (–) before to solve it formally for $m^2/T^2$; otherwise one loses the latter AF-matching properties.</footnote>
Eq. (2.26) gives almost identical results as in Fig. 5. Note also that in both Eqs. (5.1) and (5.2) the running $g(M)$ exactly cancels the M-dependence at $O(g^2)$, as easily checked using Eqs. (2.25), (2.13), and (A2).

As seen in Fig. 5 the RGOPT pressure with the (MOP or RG) $\bar{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.<fnref n="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 $\mu \neq 0$, that we omit to illustrate. We will see below that the more elaborate untruncated RGOPT pressure, accounting for higher orders in $\bar{m}(g)$, has a more drastically improved scale dependence, which is a main expected RGOPT feature.

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

### 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(\bar{m}, g, T, \mu)$ obtained from our first $\bar{m}_{MOP}$ prescription, defined by solving Eqs. (4.2), (4.9) (as explained in details in Sec. 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 $\pi T \le M \le 4\pi 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 \sim 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

<figure><img src="image_7.png" /><figcaption>FIG. 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 \le M \le 4\pi T$ at $\mu_B = 0$.</figcaption></figure>

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

### 2. Alternative RG prescription

Similarly to Fig. 6, we illustrate in Fig. 7 the exact NLO RGOPT pressure as obtained from the alternative $\bar{m}_{RG}$ prescription defined from solving Eqs. (4.10) and (4.13) (explained in details in Sec. IV D). As is seen the RGOPT reduction of remnant scale dependence is even more substantial than for the previous $\bar{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. Despite the numerically quite different MOP and RG dressed mass

<figure><img src="image_8.png" /><figcaption>FIG. 7. Same captions as for Fig. 6 but with the RGOPT pressure obtained from alternative $\bar{m}_{RG}$ prescription Eqs. (4.10), (4.13).</figcaption></figure>

<footnote n="11">¹¹We 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. [31], specially at very low $T$. This is due to considering here only its pure quark contributions, and partly also from using the exact Eq. (2.27) instead of a more approximate two-loop running expression used in [31].</footnote>
<figure><img src="image_9.png" /><figcaption>FIG. 8. $\Delta P/P \equiv P(M = 4\pi T)/P(M = \pi T) - 1$ as a function of temperature (for $\mu_B = 0$) for the different NLO RGOPT prescriptions compared to standard NLO pQCD, with scale dependence $\pi T \le M \le 4\pi T$.</figcaption></figure>

(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 and 7.

### C. Hot and dense quark matter

We now consider 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, \mu \neq 0$ case. As a representative physical value we illustrate our results for

<figure><img src="image_10.png" /><figcaption>FIG. 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} \le M \le 4\pi\sqrt{T^2 + \mu^2/\pi^2}$ at $\mu_B = 1.2$ GeV.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 10. Same captions as in Fig. 9 with alternative NLO RG prescription.</figcaption></figure>

$\mu_B = 1.2$ GeV. For the renormalization scale variation range we take as is common $\pi\sqrt{T^2 + \mu^2/\pi^2} \le M \le 4\pi\sqrt{T^2 + \mu^2/\pi^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 Figs. 9 and 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 renormal 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.

### 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 (nonlocal) 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_D^4\alpha_S$ contributions, rather than expanded in $m_D^2/T^2$ up to order $\alpha_S^{5/2}$, 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],
<figure><img src="image_12.png" /><figcaption>FIG. 11. RGOPT $P(\bar{m}(g^2))$ plus NLO $P_g^{PT}$ pressure as function of $T$ (green band) compared to $(N^3LO, g^3 \ln g)$ pQCD (light blue band), with scale dependence $\pi T \le M \le 4\pi T$, and to lattice data [50–52] at $\mu_B = 0$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 12. Full NLO RGOPT (MOP prescription) plus NLO $P_g^{PT}$ pressure as function of $T$ (grey band) compared to $(N^3LOg^3 \ln g)$ pQCD (light blue band), with scale dependence $\pi T \le M \le 4\pi T$, and to lattice data [50–52] at $\mu_B = 0$.</figcaption></figure>

$$ \frac{P_g^{PT}}{P_{g,SB}} = 1 - \frac{15}{4} \left( \frac{g}{4\pi^2} \right) + \mathcal{O}(g^2), \qquad (5.3) $$

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. (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} \rightarrow P_{g,SB} + P_{g,SB}$.<fnref n="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 reexpanded approximation for $\bar{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^3LO$ pQCD, which expression is taken from Ref. [16], and to available LQCD results from Refs. [50–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^3LO$ pQCD [16], and to LQCD

results [50–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 $\mathcal{O}(g^3 \ln g)$ pQCD[16] 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\pi T$, for temperatures as low as $T \sim 0.25$ GeV up to $T = 1$ GeV, the highest value considered in [50]. (More precisely let us mention that for the five

<figure><img src="image_14.png" /><figcaption>FIG. 13. Full NLO RGOPT (RG prescription) plus NLO $P_g^{PT}$ pressure (brown band) compared to $N^3LOg^3 \ln g$ pQCD (light blue band), NLO HTLpt (light green band) and NNLO HTLpt (light red band), with scale dependence $\pi T \le M \le 4\pi T$, and to lattice data [50–52] at $\mu_B = 0$.</figcaption></figure>

<footnote n="12">¹²As 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. (2.9), while in Figs. 11–13 below $P_{SB} \equiv P_{g,SB} + P_{g,SB}$.</footnote>
<figure><img src="image_15.png" /><figcaption>FIG. 14. NLO RGOPT (RG prescription) trace anomaly $\Delta \equiv \epsilon - 3P$ (including $\Delta_g^{PT}$) (brown band) compared to lattice data [50–52]. 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).</figcaption></figure>

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 \text{ GeV} \lesssim T \lesssim 1 \text{ 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.<fnref n="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 Sec. 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 \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 gm^4 s_2$ subtraction term of Eq. (2.16), that is formally of NLO $O(g)$.<fnref n="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

<figure><img src="image_16.png" /><figcaption>FIG. 15. Full NLO RGOPT (MOP prescription) plus NLO $P_g^{PT}$ pressure (grey band), as function of the quark chemical potential $\mu$ for $T = 0.3$ GeV, compared to $N^3LOg^3 \ln g$ pQCD (light blue band) and NNLO HTLpt (light red band), with scale dependence $\pi(T^2 + \mu^2/\pi^2)^{1/2} \le M \le 4\pi(T^2 + \mu^2/\pi^2)^{1/2}$.</figcaption></figure>

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,

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

(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 [50–52] only. A very good agreement with LQCD results of [50,51] is obtained for $0.3 \text{ GeV} \lesssim T \lesssim 1 \text{ 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 \lesssim 0.25 \text{ GeV}$.

As a last fairly different illustration, we show the NLO RGOPT pressure as function of the quark chemical potential $\mu$, for our two MOP and RG prescriptions respectively in Figs. 15 and 16, for a fixed relatively low temperature

<footnote n="13">¹³We show LQCD data as given in publicly available files [50–52] that do not include systematic uncertainties.</footnote>

<footnote n="14">¹⁴This variant is the next order analogue of including the NLO coefficient $s_1 \neq 0$ within LO RGOPT; see e.g., Eq. (3.3).</footnote>
<figure><img src="image_17.png" /><figcaption>FIG. 16. Same captions as Fig. 15 but for $\bar{m}_{RG}$ prescription (brown band).</figcaption></figure>

$T = 0.3$ GeV, thus definitely above the ($\mu = 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 $\mu$ values, and exhibits a more regular behavior at low $\mu$. 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.

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:

(i) 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 \sim \#gT^2$, and HTLpt is also based on high temperature expansions, the latter uncanceled 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].

(ii) 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.

(iii) Finally, remark that upon choosing a variational mass prescription Eq. (2.21) in HTLpt (as was done e.g., in [31,32]), nonreal $\bar{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 $\bar{m}_D$ arbitrary mass by a perturbative thermal mass [30,33], and taking the quark mass $\bar{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 $\pi T \le M \le 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, $O(g^2)$.

Our two possible 1'MOP" or "RG" prescriptions reflect the often nonuniqueness 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\pi T$ results
are in remarkable agreement with the LQCD results [50,51] for temperatures as low as $T \gtrsim 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 \gtrsim 0.25$ GeV ~ $1.5T_{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 $\lesssim T \lesssim 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,<fnref n="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.

## ACKNOWLEDGMENTS

We thank Peter Petreczky for bringing the results of Ref. [52] 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 various 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., [14,44]) of the basic thermal integrals defined in Eqs. (2.2)–(2.4),

$$2J_1(T \gg m, \mu = 0) \\ \approx \frac{7\pi^2}{180} - \frac{m^2}{12T^2} + \frac{m^4}{T^4} \left[ \frac{3}{4} - \ln\left(\frac{me^{\gamma_E}}{\pi T}\right) \right] + \mathcal{O}\left(\frac{m^6}{T^6}\right), \quad (A1)$$

$$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] + \mathcal{O}\left(\frac{m^4}{T^4}\right), \quad (A2)$$

and the more complicated genuine two-loop integral $J_3$ of Eq. (2.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) + \mathcal{O}\left(\frac{m^2}{T^2}\right) \\ \approx -0.00129532 + \mathcal{O}\left(\frac{m^2}{T^2}\right), \quad (A3)$$

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 $\bar{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 C, and IV D. Note that using directly the MOP Eq. (4.2) or the RG Eq. (4.10) 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.2) 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$. Once $M, g(M)$ are taken to be $T, \mu$-dependent via the perturbative running coupling Eq. (2.27), Eq. (4.2) becomes a function of $m/T$ and $g(T/\Lambda_{\overline{M}})$. Despite the nonlinear dependence in $m/T$, at the level of Eq. (4.2) both the AF and non-AF

<footnote n="15">¹⁵The LQCD simulations in Refs. [50–52] 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. (5.4).</footnote>
<figure><img src="image_18.png" /><figcaption>FIG. 17. AF and non-AF roots of MOP Eqs. (4.2), (4.9) for $M = 2\pi T$ and for two representative $T$ values, $T = 0.5$ GeV (dashed), $T = 1$ GeV (thick).</figcaption></figure>

solutions happen to be unique in their respective existence range. This is illustrated in Fig. 17 (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 \le M \le 4\pi T)$, $-\ln(m^2/M^2) + B_{\text{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 $\bar{m}(\text{AF}) < m_0 < \bar{m}(\text{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\pi T$, one finds $m_0 \approx 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 $\mu = 0$, following the AF-matching $\bar{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 \sim \pi T$ one hits on $M \sim \Lambda_{\overline{MS}}$ where the perturbative

coupling diverges. For sizable $\mu \ne 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. (4.10), both NLO solutions are already AF-matching, giving thus a unique solution upon using the prescription Eq. (4.13). Numerically the exact $\bar{m}_{\text{RG}}$ solution of Eq. (4.10) is somewhat larger than $\bar{m}_{\text{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[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 $\propto 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 Ref. [33] (see also Ref. [32]).

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

[2] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Szabo, Nature (London) **443**, 675 (2006); Y. Aoki, S. Borsanyi, S. Durr, Z. Fodor, S. D. Katz, S. Krieg, and K. K. Szabo, J. High Energy Phys. 06 (2009) 088; S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabó (Wuppertal-Budapest Collaboration), J. High Energy Phys. 09 (2010) 073; A. Bazavov, T. Bhattacharya, M. Cheng, C. DeTar, H. T. Ding, S. Gottlieb, R. Gupta, P. Hegde, U. M. Heller, F. Karsch et al., Phys. Rev. D **85**, 054503 (2012).

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

[4] Y. Nambu and G. Jona-Lasinio, Phys. Rev. **122**, 345 (1961).

[5] M. Buballa, Phys. Rep. **407**, 205 (2005).

[6] P. Costa, M. C. Ruivo, and C. A. de Sousa, Phys. Rev. D **77**, 096001 (2008).

[7] K. Fukushima, Phys. Lett. B **591**, 277 (2004).

[8] P. Costa, M. C. Ruivo, C. A. de Sousa, and H. Hansen, Symmetry **2**, 1338 (2010).

[9] C. D. Roberts and S. M. Schmidt, Prog. Part. Nucl. Phys. **45**, S1 (2000).

[10] C. S. Fischer, Prog. Part. Nucl. Phys. **105**, 1 (2019).
[11] W. j. Fu, J. M. Pawlowski, and F. Rennecke, Phys. Rev. D **101**, 054032 (2020); F. Gao and J. M. Pawlowski, Phys. Rev. D **102**, 034027 (2020).

[12] J. Maelger, U. Reinosa, and J. Serreau, Phys. Rev. D **97**, 074027 (2018).

[13] J. P. Blaizot, E. Iancu, and A. Rebhan, in *Quark Gluon Plasma*, edited by R. C. Hwa et al. (2004), pp. 60–122 [arXiv:hep-ph/0303185]; U. Kraemmer and A. Rebhan, Rep. Prog. Phys. **67**, 351 (2004).

[14] M. Laine and A. Vuorinen, Lect. Notes Phys. **925**, 1 (2016).

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

[16] K. Kajantie, M. Laine, K. Rummukainen, and Y. Schroder, Phys. Rev. D **67**, 105008 (2003).

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

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

[19] 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); Ann. Phys. (Amsterdam) **317**, 281 (2005).

[20] 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).

[21] A. D. Linde, Phys. Lett. **96B**, 289 (1980).

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

[23] A. Okopinska, Phys. Rev. D **35**, 1835 (1987).

[24] H. Yamada, Z. Phys. C **59**, 67 (1993).

[25] A. Duncan and M. Moshe, Phys. Lett. B **215**, 352 (1988); H. F. Jones and M. Moshe, Phys. Lett. B **234**, 492 (1990).

[26] 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).

[27] J. O. Andersen and L. Kyllingstad, Phys. Rev. D **78**, 076008 (2008).

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

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

[30] J. O. Andersen, M. Strickland, and N. Su, Phys. Rev. Lett. **104**, 122003 (2010); J. High Energy Phys. 08 (2010) 113.

[31] N. Haque, M. G. Mustafa, and M. Strickland, Phys. Rev. D **87**, 105007 (2013).

[32] J. O. Andersen, L. E. Leganger, M. Strickland, and N. Su, J. High Energy Phys. 08 (2011) 053; S. Mogliacci, J. O. Andersen, M. Strickland, N. Su, and A. Vuorinen, J. High Energy Phys. 12 (2013) 055; N. Haque, J. O. Andersen, M. G. Mustafa, M. Strickland, and N. Su, Phys. Rev. D **89**, 061701 (2014).

[33] N. Haque, A. Bandyopadhyay, J. O. Andersen, M. G. Mustafa, M. Strickland, and N. Su, J. High Energy Phys. 05 (2014) 027.

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

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

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

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

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

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

[40] J. L. Kneur and A. Neveu, Phys. Rev. D **101**, 074009 (2020).

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

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

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

[44] J. I. Kapusta and C. Gale, *Finite-Temperature Field Theory: Principles and Applications* (Cambridge University Press, Cambridge, England, 2006).

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

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

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

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

[49] E. V. Shuryak, Sov. Phys. JETP **47**, 212 (1978); S. A. Chin, Phys. Lett. **78B**, 552 (1978).

[50] S. Borsanyi, G. Endrodi, Z. Fodor, A. Jakovac, S. D. Katz, S. Krieg, C. Ratti, and K. K. Szabo, J. High Energy Phys. 11 (2010) 077.

[51] S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg, and K. K. Szabo, Phys. Lett. B **730**, 99 (2014).

[52] A. Bazavov, P. Petreczky, and J. H. Weber, Phys. Rev. D **97**, 014510 (2018).