# Renormalization group optimized $\lambda\phi^4$ pressure at next-to-next-to-leading order

Loïc Fernandez, Jean-Loïc Kneur

► To cite this version:

Loic Fernandez, Jean-Loïc Kneur. Renormalization group optimized $\lambda\phi^4$ pressure at next-to-next-to-leading order. Physical Review D, 2021, 104 (9), pp.096012. 10.1103/PhysRevD.104.096012. hal-03319381

HAL Id: hal-03319381

https://hal.science/hal-03319381

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 optimized $\lambda\phi^4$ pressure at next-to-next-to-leading order

Loïc Fernandez and Jean-Loïc Kneur

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

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

(Received 2 August 2021; accepted 13 October 2021; published 17 November 2021)

We investigate the renormalization group optimized perturbation theory (RGOPT) at the next-to-next-to-leading order (NNLO) for the thermal scalar field theory. From comparing three thus available successive RGOPT orders, we illustrate the efficient resummation and very good apparent convergence properties of the method. In particular, the remnant renormalization scale dependence of thermodynamical quantities is drastically improved as compared to both standard perturbative expansions and other related resummation methods, such as the screened perturbation theory. Our present results thus constitute a useful first NNLO illustration in view of NNLO applications of this approach to the more involved thermal QCD.

DOI: 10.1103/PhysRevD.104.096012

# I. INTRODUCTION

For thermodynamical quantities at equilibrium and a weakly coupled theory, one could hope at first that a perturbative expansion would give reliable results. In contrast, as is well-known, for a massless theory, infrared divergences spoil a naive perturbation theory (PT) approach in thermal field theories. Even though those infrared divergences can be efficiently resummed (see, e.g., [1–3] for reviews), it leads to nonanalytical terms in the coupling, which happen to give poorly convergent successive orders, even when pushed to the highest perturbative order available and showing increasingly sizeable remnant scale dependence. This situation is notoriously illustrated for asymptotically free thermal QCD, where lattice simulations (LS) offer a powerful genuinely nonperturbative alternative to bypassing those issues. At least so far, LS has been very successful in the description of the nonperturbative physics of the QCD phase transitions at finite temperatures and near vanishing or small baryonic densities [4]. Nevertheless, the famous numerical sign problem [5] at finite density (equivalently at finite chemical potential) prevents LS to successfully describe compressed baryonic matter at sufficiently high densities and to explore a large part of the QCD phase diagram. More generally, despite the very success of LS, it is still highly desirable to explore more analytical improvements of thermal PT. Accordingly, many efforts have been devoted in the past to overcome the generically observed issues of poor PT convergence. Apart from the

most important case of QCD, the above mentioned behavior is generic for any thermal quantum field theory, and typically, the scalar $\lambda\phi^4$ interaction is often used as a simpler model to study new alternative approaches. Although the $\lambda\phi^4$ model is not asymptotically free, it shares some important features with thermal QCD. Typically, the dynamical generation of a thermal screening mass $m_D \sim \sqrt{\lambda}T$ impacts the relevant expansion of thermodynamical quantities, such as the pressure that exhibit weak expansion terms $\lambda^{(2p+1)/2}$, $p \ge 1$.

Various approximations attempting to more efficiently resum thermal perturbative expansions have been developed and refined over the years, typically the screened perturbation theory (SPT) [6,7], the nonperturbative renormalization group (NPRG) [8] approach, the two-particle irreducible (2PI) formalism [9–11], or other approaches [12]. In particular, for the $\lambda\phi^4$ model, SPT essentially redefines the weak expansion about a quasiparticle mass, avoiding in this way infrared divergences from the start. It has been investigated up to three-loop [13–15] and even four-loop orders [16]. SPT may also be viewed as a particular case, in the thermal context, of the so-called optimized perturbation theory (OPT),<fnref n="1" /> in which more generally the (thermal or nonthermal) perturbative expansion is redefined about an unphysical test mass parameter, fixed by a variational prescription, that provides a resumination of perturbative expansion. The generalization of SPT tailored to treat the much more involved thermal gauge theories [20], the hard thermal loop perturbation theory (HTLpt) [21], has been pushed to three-loop order [22,23]. The three-loop results [23] agreement with LS is quite

<footnote>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>

<footnote n="1">¹OPT and its many variants appear under different names in the literature [17–19].</footnote>
remarkable down to about twice the critical temperature, for a renormalization scale choice ~2πT. However, both SPT and HTLpt exhibit a very sizeable remnant scale dependence at NNLO, that definitely call for further improvement.

More recently, OPT at vanishing temperatures and densities was extended to the so-called renormalization group optimized perturbation theory (RGOPT) [24,25]. The basic novelty is that it restores perturbative RG invariance at all stages of calculations, in particular when fixing the variational mass parameter, by solving the (mass) optimization prescription consistently with the RG equation. At vanishing temperatures and densities, it has given precise first principle determinations [25] of the basic QCD scale (Λ<sub>MS</sub>) or related coupling α<sub>S</sub> and of the quark condensate [26,27]. The RGOPT was extended at finite temperatures for the scalar λφ<sup>4</sup> in [28,29] and for the nonlinear sigma model (NLSM) in [30], showing how it substantially reduces the generic scale dependence and convergence problem of thermal perturbation theories at increasing perturbative orders. More recently, the RGOPT in the quark sector contribution to the QCD pressure was investigated at NLO, for finite densities and vanishing temperatures [31], and at finite temperature and density [32,33], leading to drastic improvements with respect to both perturbative QCD and HTLpt, specially at nonzero temperature. In the present work, as a first step to investigate the RGOPT in thermal theories beyond NLO, we explore the three-loop order (NNLO) for the technically simpler scalar λφ<sup>4</sup> model pressure. We investigate the stability and convergence properties of the method, also assessing the remnant renormalization scale dependence improvement as compared to standard PT, and to SPT.

The paper is organized as follows. In Sec. II, we briefly review the standard thermal perturbative pressure of the λφ⁴ model up to NNLO, to set our conventions and basic expressions that serve as a starting point for our construction. In Sec. III, we recall the main ingredients of the RGOPT construction and some previous NLO results from [28,29]. Our main new NNLO results are derived and illustrated in Sec. IV, while finally some conclusions and outlook are given in Sec. V. Some additional technical ingredients can be found in three appendixes.

<figure><img src="image_3.png" /><figcaption>FIG. 1. Free energy diagrams up to NNLO in λφ⁴ model.</figcaption></figure>

## II. REVIEW OF STANDARD (MASSIVE) THERMAL PERTURBATIVE EXPANSION

We consider the Lagrangian for one neutral scalar field with a quartic interaction,

$$ \mathcal{L} = \frac{1}{2} \partial_{\mu} \phi \partial^{\mu} \phi - \frac{m^2}{2} \phi^2 - \frac{\lambda}{4!} \phi^4, \quad (2.1) $$

where a generic mass term m is arbitrary at this stage.

### A. Free energy up to NNLO

We first recall the result [2,13] for the two-loop free energy (equivalently, minus the pressure) including a mass term, corresponding to the first two graphs<fnref n="2" /> in Fig. 1,

$$ \mathcal{F}_0 = \frac{1}{2} \int_P^c \ln(\mathbf{P}^2 + m_0^2) + \frac{\lambda_0}{8} \left( \int_P^c \frac{1}{\mathbf{P}^2 + m_0^2} \right)^2 + \mathcal{F}_0^{2l,ct}, \quad (2.2) $$

where in the imaginary time formalism $\mathbf{P}^2 = \omega_n^2 + \mathbf{p}^2$, $\omega_n = 2\pi Tn$ ($n = 0, \pm 1, ...$) represents the bosonic Matsubara frequencies. The sum integral in Eq. (2.2) is defined as usual as the sum over Matsubara frequencies times remaining integration over the three momentum, using dimensional regularization and the $\overline{\text{MS}}$ renormalization scheme,

$$ \int_P^c \equiv T \sum_n \left( \frac{\mu^2 e^{\gamma_E}}{4\pi} \right)^\epsilon \int \frac{d^{3-2\epsilon} \mathbf{p}}{(2\pi)^{3-2\epsilon}}. \quad (2.3) $$

The three-loop contributions involves the basic last two graphs in Fig. 1 and read

$$ \mathcal{F}_0^{3l} = -\frac{\lambda_0^2}{48} \left[ 3 \left( \int_P^c \frac{1}{\mathbf{P}^2 + m_0^2} \right)^2 \int_q^c \frac{1}{(\mathbf{Q}^2 + m_0^2)^2} + \int_{pqr}^c \frac{1}{(\mathbf{P}^2 + m_0^2)(\mathbf{Q}^2 + m_0^2)(\mathbf{R}^2 + m_0^2)((\mathbf{P} + \mathbf{Q} + \mathbf{R})^2 + m_0^2)} \right] + \mathcal{F}_0^{3l,ct}. \quad (2.4) $$

We recall that the renormalization is most easily performed as follows. First, applying multiplicative renormalization to the (bare) coupling and mass in expressions above, as

$$ \lambda_0 = \lambda Z_\lambda, \quad m_0 = m Z_m, \quad (2.5) $$

where $Z_λ$, $Z_m$ are the standard coupling and mass counterparts for the massive λφ⁴ model, given for completeness, respectively, in Eqs. (B5), (B6) in Appendix B. Those are

<footnote n="2">²In Fig. 1, counterterm graphs are omitted for simplicity.</footnote>
expanded for the relevant three-loop order case here to the perturbative order $\lambda^2$. The remaining divergences are then removed by an additive renormalization from the vacuum energy counterterms [13], formally represented as $\mathcal{F}_0^{2l,ct}$, $\mathcal{F}_0^{3l,ct}$ in the above expressions, and also given for completeness in Eq. (B7) in Appendix B.

Once the mass, coupling, and vacuum energy counterterms have been accounted to cancel the original divergences, one obtains the (MS scheme) renormalized free energy. The one- and two-loop contributions read [2,13]

$$ (4\pi)^2 \mathcal{F}^{\text{NLO}} = \mathcal{E}_0 - \frac{1}{8}m^4(3+2L) - \frac{T^4}{2}J_0\left(\frac{m}{T}\right) + \frac{1}{8}\left(\frac{\lambda}{16\pi^2}\right)\left[(L+1)m^2 - T^2 J_1\left(\frac{m}{T}\right)\right]^2, \quad (2.6) $$

where $L \equiv \ln(\mu^2/m^2)$, and we explicitly separated the thermal and nonthermal contributions. Here and in all related renormalized expressions below, $\mu$ stands for the arbitrary renormalization scale introduced by dimensional regularization in the MS scheme, Eq. (2.3), and $\lambda \equiv \lambda(\mu)$. Note carefully that $\mathcal{E}_0$ in Eq. (2.6) represents a finite vacuum energy term, to be specified below, that plays a crucial role in our approach as will be reexamined in Sec. II B.

The standard (dimensionless) thermal integrals appearing in Eq. (2.6) and below are given by

$$ J_n(x) = 4 \frac{\Gamma[1/2]}{\Gamma[5/2-n]} \int_0^\infty dt \frac{t^{4-2n}}{\sqrt{t^2 + x^2} e^{\sqrt{t^2+x^2}} - 1}, \quad (2.7) $$

where $t = p/T$ and $x = m/T$. Different integrals can be easily related by employing derivatives, such as

$$ J_{n+1}(x) = -\frac{1}{2x} \frac{\partial}{\partial x} J_n(x). \quad (2.8) $$

Also, a high-T expansion [2], such as

$$ J_0(x) \approx \frac{16}{45}\pi^4 - 4\frac{\pi^2}{3}x^2 + 8\frac{\pi}{3}x^3 + x^4 \left[ \ln\left(\frac{x}{4\pi}\right) + \gamma_E - \frac{3}{4} \right] + \mathcal{O}(x^6), \quad (2.9) $$

is often a rather good approximation as long as $x \lesssim 1$, i.e., $m \ll T$.

Next, the NNLO contribution involves the additional genuine massive three-loop second integral in Eq. (2.4), first calculated in [34]. After algebra, the complete three-loop contribution can be expressed as [13]

$$ F_{3l} = -\frac{1}{48} \left(\frac{\lambda}{16\pi^2}\right)^2 \left[ \begin{aligned} & m^4 \left( 5L^3 + 17L^2 + \frac{41}{2}L - 23 - \frac{23}{12}\pi^2 + 2\zeta(3) + c_0 + 3(L+1)^2 J_2\left(\frac{m}{T}\right) \right) \\ & - m^2 T^2 J_1\left(\frac{m}{T}\right) \left( 12L^2 + 28L - 12 - \pi^2 - 4c_1 + 6(L+1)J_2\left(\frac{m}{T}\right) \right) \\ & + T^4 \left( 3\left(3L+4+J_2\left(\frac{m}{T}\right)\right) J_1^2\left(\frac{m}{T}\right) + 6K_2\left(\frac{m}{T}\right) + 4K_3\left(\frac{m}{T}\right) \right) \end{aligned} \right], \quad (2.10) $$

where

$$ c_0 = \frac{275}{12} + \frac{23}{2}\zeta(2) - 2\zeta(3) \approx 39.429, \quad c_1 = -\frac{59}{8} - \frac{3}{2}\zeta(2) \approx -9.8424, \quad (2.11) $$

and it involves two irreducible, respectively, two-loop $K_2(m/T)$ and three-loop $K_3(m/T)$ integrals, given explicitly in Ref. [34], reproduced for self-containedness in Eqs. (C1), (C4) in Appendix C.

To complete this subsection on basic thermal perturbative expressions, for an easier latter reference and comparison purpose, we also recall the expression of the (massless) PT pressure [1,35,36] up to NNLO,

$$ \frac{P}{P_0} = 1 - \frac{5}{4} \left( \frac{\lambda}{16\pi^2} \right) + 5 \frac{\sqrt{6}}{3} \left( \frac{\lambda}{16\pi^2} \right)^{3/2} + \frac{15}{4} \left( \frac{\lambda}{16\pi^2} \right)^2 \left( \ln \frac{\mu}{2\pi T} + 0.40 \right) + \mathcal{O}(\lambda^{5/2}), \quad (2.12) $$

where $P_0 = (\pi^2/90)T^4$ is the ideal bosonic gas pressure.
## B. Perturbatively RG-invariant massive free energy

At this stage, an important feature is that the (massive) free energy is lacking RG invariance. Namely, applying to Eq. (2.6), with $\mathcal{E}_0 = 0$, the standard RG operator,

$$ \mu \frac{d}{d\mu} = \mu \frac{\partial}{\partial \mu} + \beta(\lambda) \frac{\partial}{\partial \lambda} + \gamma_m(\lambda) m \frac{\partial}{\partial m}, \quad (2.13) $$

with $\beta(\lambda)$, $\gamma_m(\lambda)$ given in Eqs. (B1), (B2), yields a remnant contribution of leading (one-loop) order: $-(m^4/2) \ln \mu$, for arbitrary $m$. Indeed, the latter term is not compensated by the lowest orders terms from $\beta(\lambda)$ or $\gamma_m(\lambda)$ in Eq. (2.13), those being at least of next order $O(\lambda m^4)$. This is a manifestation that perturbative RG invariance generally occurs from cancellations between terms from the RG equation at an order $\lambda^k$ and the explicit $\mu$ dependence at the next order $\lambda^{k+1}$. Nevertheless, perturbative RG invariance can easily be restored by adding the finite vacuum energy term $\mathcal{E}_0$ to the action, although this term is usually ignored and minimally set to zero in the (thermal) literature [13,21,22]. Following Refs. [28,29], the easiest way to construct a (perturbatively) RG-invariant renormalized vacuum energy is to determine $\mathcal{E}_0$ order by order as a perturbative series from the reminder of Eq. (2.13) applied on the non-RG-invariant finite part of Eq. (2.6),

$$ \mu \frac{d}{d\mu} \mathcal{E}_0(\lambda, m) \equiv -\text{Remnant}(\lambda, m) \\ = -\mu \frac{d}{d\mu} [\mathcal{F}_0(\mathcal{E}_0 \equiv 0)|_{\text{finite}}]. \quad (2.14) $$

Accordingly, $\mathcal{E}_0$ has the form,

$$ \mathcal{E}_0(\lambda, m) = -\frac{m^4}{\lambda} \sum_{k \ge 0} s_k \lambda^k, \quad (2.15) $$

where the coefficients $s_k$ can be determined at successive orders from knowing the (single) powers of $\ln \mu$ at order $k+1$ (or equivalently, the single poles in $1/\epsilon$ of the unrenormalized expression) [25]. This procedure leaves non-RG-invariant remnant terms of higher orders, which may be treated similarly once higher order terms are considered. Explicitly, we obtain [28]

$$ \begin{aligned} s_0 &= \frac{1}{2(b_0 - 4\gamma_0)} = 8\pi^2 \\ s_1 &= \frac{(b_1 - 4\gamma_1)}{8\gamma_0(b_0 - 4\gamma_0)} = -1, \end{aligned} \quad (2.16) $$

and similarly, for the next orders presently relevant,

$$ \begin{aligned} s_2 &= \frac{96\pi^2(b_0 - 128\pi^2((1 + 4s_1)\gamma_1 - s_0(b_2 - 4\gamma_2)) - 41)}{12288\pi^4(b_0 + 4\gamma_0)} \\ &= \frac{23 + 36\zeta[3]}{480\pi^2} \approx 0.01399, \\ s_3 &= \frac{-709 + 12\pi^4 - 2628\zeta(3) - 5400\zeta(5)}{720(16\pi^2)^2}. \end{aligned} \quad (2.17) $$

The explicit RG coefficients in intermediate expressions emphasize the general form of these results, while the $s_3$ expression is specific to the $N=1$ $\lambda\phi^4$ model.

We stress that the previous construction, being only dependent on the renormalization procedure, does not depend on temperature-dependent contributions: at arbitrary perturbative orders, the $s_k$ coefficients are determined from the $T=0$ contributions only. Indeed, as Eq. (2.14) suggests, its rhs precisely defines the vacuum energy anomalous dimension that has been calculated even to the five-loop order for the general $O(N)$ scalar model [37]. Our independent results for the $s_k$ are fully consistent with [37]. A subtlety is that according to Eq. (2.15), $s_k$ is strictly required for perturbative RG invariance at the order $\lambda^k$, but contributes at the order $\lambda^{k-1}$. So at an order $\lambda^k$, one may choose minimally to include only $s_0, \dots, s_k$, or more completely include also $s_{k+1}$, incorporating in this way a higher order RG dependence within the resulting expression.

## III. RGOPT $\lambda\phi^4$ PRESSURE AT LO AND NLO

In this section, we briefly recall the RGOPT construction, as investigated for the $\lambda\phi^4$ model in [28,29] at LO and NLO, before extending our approach to the technically more involved NNLO. We also underline some important features that were not plainly discussed in [28]. After restoring in a first stage perturbative RG invariance of the massive free energy, leading to the crucial additional term in Eq. (2.15), one performs on the resulting complete expression, the variational modification, according to

$$ \mathcal{F}_0(m^2 \rightarrow m^2(1 - \delta)^{2a}, \lambda \rightarrow \delta\lambda), \quad (3.1) $$

where $m$ is from now an arbitrary variational mass, and the crucial role of the exponent $a$ will be specified just next. One then reexpands Eq. (3.1) to successive orders, $\delta^k$, at the same order as the original perturbative expression, and set $\delta \rightarrow 1$ afterwards. This leaves a remnant $m$ dependence at any order $k$, that may be conveniently fixed by a stationarity prescription [19],

$$ \frac{\partial \mathcal{F}_0^{(k)}}{\partial m} (m, \lambda, \delta = 1)|_{m=\tilde{m}} \equiv 0, \quad (3.2) $$

thus determining a dressed mass $\tilde{m}(\lambda)$ with a “nonperturbative” (all order) $\lambda$ dependence. In OPT [17] or similarly SPT [7] applications, the linear $\delta$ expansion has been
mostly used, i.e., assuming $a = 1/2$ in Eq. (3.1), and that corresponds to the “add and subtract a mass” intuitive prescription with one mass being treated as an interaction term. Yet it was pointed out that the rather drastic modification implied by Eq. (3.1) is generally not compatible with RG invariance [28]: in contrast, $a$ can be uniquely fixed [25,28] by (re)imposing RG invariance now for the variationally modified perturbative expansion. Once combined with Eq. (3.2), the RG Eq. (2.13) takes the massless form,

$$ \left[ \mu \frac{\partial}{\partial \mu} + \beta(\lambda) \frac{\partial}{\partial \lambda} \right] \mathcal{F}_0^{(k)}(m, \lambda, a, \delta = 1) = 0, \quad (3.3) $$

that at leading RG order uniquely fixes [25,28]

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

simply in terms of the universal (renormalization scheme independent) first order RG coefficients.<fnref n="3" />

At higher orders, Eq. (3.3) is no longer exactly fulfilled and can be thus used to determine an RG-compatible dressed mass, $\tilde{m}_{RG}(\lambda, T)$, as a possible alternative to the OPT Eq. (3.2). Importantly, Eq. (3.4) guarantees in addition that either Eq. (3.2) or Eq. (3.3) have at least one (essentially unique) solution matching the $T = 0$ perturbative behavior [25,28] for $\lambda \to 0$, i.e., infrared freedom in the present case: $\lambda(\mu \ll m) \approx [b_0 \ln(m/\mu)]^{-1}$.

## A. LO RGOPT

For a simple illustration, let us briefly recall the one-loop results [28,29], thus considering the LO term in Eq. (2.2), including solely the first order subtraction term, $\mathcal{E}_0 = -(m^4/\lambda)s_0$, with $s_0$ given in Eq. (2.16). Performing (3.1), expanding to leading order $\delta^0$ consistently, and taking afterwards $\delta \to 1$ yields

$$ (4\pi)^2 \mathcal{F}_0^{\delta^0} = -m^4 \left[ \frac{1}{2b_0\lambda} + \left( \frac{3}{8} + \frac{1}{4}\ln\frac{\mu^2}{m^2} \right) \right] - \frac{T^4}{2} J_0\left(\frac{m}{T}\right). \quad (3.5) $$

At this leading order, Eq. (3.3) is satisfied exactly, so that only Eq. (3.2) can determine a nontrivial dressed thermal mass. It is convenient to introduce first the one-loop renormalized self-energy, including all the relevant $T$ dependence, $\Sigma_R$, explicitly [1,13],

$$ \Sigma_R(m) = \gamma_0 \lambda \left[ m^2 \left( \ln \frac{m^2}{\mu^2} - 1 \right) + T^2 J_1 \left( \frac{m}{T} \right) \right]. \quad (3.6) $$

Then the exact solution of Eq. (3.2), using Eq. (2.8), is given by the self-consistent gap equation,

$$ \tilde{m}^2 = \frac{1}{2} \left( \frac{b_0}{\gamma_0} \right) \Sigma_R(\tilde{m}^2) = (4\pi)^2 b_0 \Sigma_R(\tilde{m}^2). \quad (3.7) $$

Accordingly, $\tilde{m}$ is exactly (one-loop) RG invariant: $\lambda \equiv \lambda(\mu)$ being given by the “exact” (one-loop) running,

$$ \frac{1}{\lambda(\mu)} = \frac{1}{\lambda(\mu_0)} - b_0 \ln \frac{\mu}{\mu_0}, \quad (3.8) $$

the free energy Eq. (3.5), and therefore, $\tilde{m}$ in Eq. (3.7), only depend on the combination $1/(b_0\lambda(\mu)) + 1/2 \ln(\mu^2/m^2)$ that is $\mu$ independent. Let us remark at this stage that some interesting qualitative similarities were observed [28] between those (one-loop) RGOPT results and the (two-loop) 2PI resummation approach in [9,10], although our construction is basically very different.

To get more insight on some properties of the solution of Eq. (3.7), one may conveniently use the high-temperature expansion of the relevant $J_n(x)$, from Eq. (2.9) with $x \equiv m/T$.<fnref n="4" /> Accordingly, Eq. (3.7) simplifies to a quadratic equation for $x$, with a unique physical ($x > 0$) solution,

$$ \begin{aligned} \tilde{x} &= \frac{\tilde{m}^{(1)}}{T} = \pi \frac{\sqrt{1 + \frac{2}{3} \left( \frac{1}{b_0 \lambda} + L_T \right) - 1}}{\frac{1}{b_0 \lambda} + L_T} \\ &\approx \pi \sqrt{\frac{2}{3} b_0 \lambda - \pi b_0 \lambda + \mathcal{O}(\lambda^{3/2})}, \end{aligned} \quad (3.9) $$

with $L_T \equiv \ln[\mu e^{\gamma_E}/(4\pi T)]$. Perturbatively at high temperature, the variationally determined mass Eq. (3.9) has the form of a screening mass, $\tilde{m}^2 \sim \mathcal{O}(\lambda T^2)$, but note that this variational mass parameter is unrelated to the physical Debye mass [36] definition. The corresponding one-loop RGOPT pressure, from Eq. (3.5), reads

$$ \frac{P^{(1)}}{P_0}(\tilde{x}) = 1 - \frac{15}{4\pi^2}\tilde{x}^2 + \frac{15}{2\pi^3}\tilde{x}^3 + \frac{45}{16\pi^4}\left(\frac{1}{b_0\lambda} + L_T\right)\tilde{x}^4 + \mathcal{O}(\tilde{x}^6). \quad (3.10) $$

Equations (3.7)–(3.10) clearly involve an all order dependence in $\lambda$, and RG invariance is once more manifest since from Eq. (3.8), $1/(b_0\lambda(\mu)) + L_T$ is $\mu$ independent. Accordingly, Eqs. (3.9) and (3.10) only depend on the

<footnote n="3">³At higher orders, one could generalize the interpolation $(1 - \delta)^a$ with $\delta^2$ and higher order terms without spoiling the crucial RG properties guaranteed from Eq. (3.4). But this would involve extra arbitrary variational parameters with no compelling reasons. We thus keep the simpler form Eq. (3.4) at successive orders for sensible comparisons.</footnote>

<footnote n="4">⁴At the one-loop order, this approximation is valid at the 0.1% level even for $x \lesssim 1$, sufficient for our purpose since the RGOPT one-loop solution $\tilde{m}/T$ happens to always lie in this range.</footnote>
<figure><img src="image_4.png" /><figcaption>FIG. 2. The graphs being resummed at first nontrivial RGOPT order.</figcaption></figure>

single parameter $b_0\lambda(\mu_0)$, where $\mu_0$ is some reference scale, typically $\mu_0 = 2\pi T_0$.

Expanding perturbatively Eq. (3.10), one obtains for the first few orders,

$$ \frac{P^{(1)}}{P_0} = 1 - \frac{5}{4}\alpha + 5\frac{\sqrt{6}}{3}\alpha^{3/2} + \frac{5}{4}(L_T - 6)\alpha^2 + \mathcal{O}(\alpha^{5/2}), \qquad (3.11) $$

where $\alpha \equiv b_0\lambda$. Note in particular that Eq. (3.11) contains the nonanalytic term $\sim \lambda^{3/2}$, that originates from the bosonic zero mode resummation, but here is readily obtained from RG properties. Expanding Eq. (3.9) at higher orders, it is easily seen that it entails the nonanalytic terms $\lambda^{(2p+1)/2}$, $p \ge 1$ at all orders. Incidentally, it is worth mentioning that the higher orders beyond Eq. (3.11) obtained from Eq. (3.9) correctly reproduce all orders of the $O(N)$ scalar model large $N$ results [e.g., Eq. (5.8) of [38]] [once including higher order $O(x^6)$ terms, not given in Eq. (2.9)], as can be checked upon identifying the correct large-$N$ $b_0 = 1/(16\pi^2)$ value [38]. Accordingly, although the LO RGOPT is essentially built on the very first one-loop graph of Fig. 1 augmented by the optimized RG construction as above described, it happens to correctly resum the whole set of “foam” graphs as illustrated in Fig. 2.

## B. NLO RGOPT

The NLO (two-loop) $O(\delta^1)$ contribution to the free energy, for $\delta = 1$, takes a rather compact form in terms of $\Sigma_R$ in Eq. (3.6),

$$ \mathcal{F}_0^{\delta^1} = \frac{\mathcal{E}_0^{\delta^1}}{(4\pi)^2} + \frac{T}{2} \int_{\mathbf{p}} \ln(\omega_n^2 + \omega_{\mathbf{p}}^2)|_{\text{finite}} \\ - \left(\frac{2\gamma_0}{b_0}\right) \frac{m^2}{\lambda} \Sigma_R + \frac{\Sigma_R^2}{2\lambda}, \qquad (3.12) $$

where the subscript on the integral term means to taking its corresponding finite (renormalized) expression. Equation (2.15) gives explicitly,

$$ \mathcal{E}_0^{\delta^1} = -m^4 \left( \frac{1}{3b_0\lambda} + \frac{s_1}{3} \right). \qquad (3.13) $$

The exact two-loop OPT and RG Eqs. (3.2) and (3.3) can be written compactly as [28]$^5$

$$ f_{\text{OPT}}^{\text{NLO}} = \frac{2}{3}h\left(-s_1 - \frac{1}{b_0\lambda}\right) + \frac{2}{3}S + \frac{\Sigma'_R}{2m}\left(S - \frac{1}{3\lambda}\right) = 0, \qquad (3.14) $$

$$ f_{\text{RG}}^{\text{NLO}} = h \left[ \frac{1}{6} + \left( \frac{b_1}{3b_0} - S \right) \lambda \right] + \frac{1}{2} \beta^{(2)}(\lambda) S^2 = 0, \qquad (3.15) $$

with $h \equiv (4\pi)^{-2}$, $\beta^{(2)}(\lambda) = b_0\lambda^2 + b_1\lambda^3$, and we introduced the reduced (dimensionless) self-energy,

$$ S(m, \mu, T) \equiv \Sigma_R / (m^2 \lambda), \qquad (3.16) $$

with, from Eq. (3.6),

$$ \begin{aligned} \Sigma'_R &\equiv \partial_{m^2}(\Sigma_R) = \lambda(S + m^2 S') \\ &= \gamma_0 \lambda[\ln(m^2/\mu^2) - J_2(m/T)]. \end{aligned} \qquad (3.17) $$

At this stage, at NLO, in principle, we could use three different possible prescriptions to obtain a thermally dressed mass $\tilde{m}(\lambda, T)$ as a function of the coupling, as investigated in details in [28]: either the OPT Eq. (3.14) or alternatively, the (massless) RG Eq. (3.15) or else the full RG, Eq. (2.13). The latter is not an independent equation, being a linear combination of Eqs. (3.14) and (3.15),

$$ f_{\text{full RG}}^{\text{NLO}} \equiv f_{\text{RG}} + 2\gamma_m(\lambda)f_{\text{OPT}} = 0. \qquad (3.18) $$

We speculate that if one could calculate to all orders, the solutions of those different prescriptions would presumably converge towards a unique, nonperturbatively dressed mass $\tilde{m}(g, T)$ [as it happens [24] in the large-$N$ limit of the $O(N)$ Gross-Neveu model, where the original perturbative series is known to all orders]. But due to the inherent perturbative truncations, those prescriptions give formally different solutions, thus providing useful variants of the method. The resulting NLO solutions for $\tilde{m}/T$ and $P/P_0$, once reexpanded, are perturbatively consistent with Eqs. (3.9), (3.10) for the first order term but contain modifications at higher orders [28]. More precisely, in the standard OPT/SPT prescription, $a = 1/2$ in Eq. (3.1) typically, the solution of Eq. (3.2) would be a priori very different from the RG ones of Eq. (3.3) or Eq. (3.18). But a remarkable consequence of Eq. (3.4) is that all the three solutions have the same leading term (in $\lambda$) given by the LO term in Eq. (3.9). However, for rather large $g(2\pi T) \lesssim 1$ reference coupling and $\mu = 4\pi T$,

<footnote>$^5$Note a typo in Eq. (7.2) of Ref. [28], corrected in Eq. (3.14), that does not affect any numerical results.</footnote>
<figure><img src="image_5.png" /><figcaption>FIG. 3. LO (dotted lines) and NLO (dashed, green) RGOPT pressure $P/P_0(g \equiv \sqrt{\lambda/24})$ versus $O(g^2)$ and $O(g^{3/2})$ PT, and NLO SPT (for two different prescriptions, see main text). The different bands give the scale dependence for $\pi T \le \mu \le 4\pi T$.</figcaption></figure>

the OPT Eq. (3.14) no longer gives a real solution: in fact, the possible occurrence of nonreal solutions at higher orders is a recurrent burden of such a variational approach, quite generically expected from the nonlinear $m$ dependence if requiring to solve $\tilde{m}(\lambda, T)$ exactly. Alternatively, using Eq. (3.15) to determine $\tilde{m}(\lambda, T)$ gives an unphysical solution at NLO [28], being driven towards the NLO UV fixed point at $\lambda = -b_0/b_1$ (an artifact of the scalar model two-loop beta function approximation due to $b_1 < 0$). These features lead us to rather consider the third option at NLO, taking the full RG Eq. (2.13), explicitly Eq. (3.18) at NLO, that happens to give real solutions, at least in a much larger range of coupling values.<fnref n="6" />

The NLO pressure $P/P_0$ with exact $T$ dependence obtained from Eq. (3.18), as a function of the reference coupling $g \equiv \sqrt{\lambda(2\pi T)/24}$, is illustrated in Fig. 3, with a scale dependence $\pi T \le \mu \le 4\pi T$ using the exact two-loop running Eq. (B8). It is compared with the LO RGOPT, Eq. (3.10), and importantly, also with the PT pressure, Eq. (2.12), respectively, truncated at the lowest $O(\lambda) \sim g^2$ and next $O(\lambda^{3/2}) \sim g^3$ orders. We recall that the latter nonanalytic terms, obtained from resummation of a certain class of individually infrared divergent (massless) graphs, are largely responsible for the poorly convergent, oscillating behavior of successive perturbative orders, as illustrated. In contrast, going from LO to NLO RGOPT appears very stable, despite the fact that both approximations also incorporate $O(\lambda^{3/2}) \sim g^3$ contributions (as well as arbitrary higher order contributions as above explained). Moreover, the reduction of the remnant scale dependence as compared to standard PT pressure is also sizeable, although a moderate residual scale dependence appears at NLO,

visible on the figure for $g \gtrsim 0.6$. A more accurate analysis [28] shows that at NLO, the remnant scale dependence reappears first at the perturbative order $\lambda^3$. As mentioned above, the exact scale invariance obtained at (one-loop) LO RGOPT (as illustrated in Fig. 3, where the dotted line has no visible width) results from the peculiar form of the exact running coupling Eq. (3.8) perfectly matching Eqs. (3.9), (3.10), since the latter only depends on $b_0$. In contrast, at NLO and beyond, $P(\tilde{m}(\lambda), T, \mu)$ inevitably has a remnant scale dependence, basically because the subtractions in Eq. (2.15) only guarantee RG invariance up to remnant higher order terms, $\sim m^4\lambda^2$ at NLO. It is a nontrivial consequence of the subsequent variational modification in Eq. (3.1) that it also preserves this RG invariance at the same level, i.e., up to basically neglected terms of formally higher orders.

The LO and NLO RGOPT pressures in Fig. 3 trivially reproduce the Stefan-Boltzmann limit for $\lambda \to 0$; see Eq. (3.10). But one can notice that beyond moderate coupling values, they differ rather importantly from the PT pressure: this is more pronounced in such a plot, often conventionally used [1] to illustrate different approximations to $P(\lambda)$, but where the $\phi^4$ model is not fully specified by fixing a physical input scale $\mu = T_0$ and corresponding $\lambda(T_0)$ value. Note however that upon expressing our results in terms of a more physical mass scale, namely solving, e.g., Eq. (3.2) for $\tilde{\lambda}(m)$ with a resulting $P(m/T)$, and inserting the physical Debye screening mass $m(\lambda)$ [36], correctly reproduces [28] the first two terms of the standard PT pressure Eq. (2.12). Yet the RGOPT pressure crucially differs from PT at higher orders; otherwise, it would merely reproduce the same PT behavior and issues. Thus, it is important to use the exact (all order) $\tilde{m}(\lambda, T)$ solution from Eqs. (3.14) or (3.15) that corresponds to the results shown in Fig. 3. (In contrast, using low orders perturbative reexpansions of the RGOPT $\tilde{m}(\lambda)$ solution leads to a behavior more similar to the PT pressure, showing large differences between successive orders and a larger scale dependence.)

In Fig. 3, we also compare with the NLO SPT results elaborated in Ref. [13]. Note that all the relevant SPT expressions can be obtained consistently by 1) discarding the vacuum energy subtraction $\mathcal{E}_0$ in Eq. (2.6); 2) taking $a = 1/2$ in Eq. (3.1), expanding the result to order $\delta$, and setting $\delta \to 1$; and finally, 3) calculating the variational mass gap, Eq. (3.2) that gives explicitly,<fnref n="7" />

$$m_{\text{SPT}}^2 = \Sigma_R, \qquad (3.19)$$

see Eq. (3.6), to be solved self-consistently for $\tilde{m}_{\text{SPT}}$. Alternatively, a simpler prescription was also used in Ref. [13], taking instead the (NLO) perturbative Debye screening mass [36],

<footnote n="6">⁶In RGOPT applications to thermal QCD [31,33], at NLO, nonreal solutions occur at increasing QCD coupling values for all the prescriptions, an issue that can be circumvented at the price of more elaborate prescriptions, based on renormalization scheme changes.</footnote>

<footnote n="7">⁷Equation (3.19) is what is called the tadpole prescription in Ref. [13].</footnote>
$$m_D^2 = \frac{\lambda}{24} T^2 \left(1 - \frac{3}{\pi} \left(\frac{\lambda}{24}\right)^{1/2}\right); \quad (3.20)$$

therefore, we illustrate these two SPT prescriptions in Fig. 3. As compared to the two lowest orders of standard PT shown, the SPT pressure is significantly more stable and with a better remnant scale dependence that reflects its more elaborate resummation properties. The SPT pressure values obtained from the two prescriptions are quite close, but using the variational mass gap gives a much better remnant scale dependence than using the PT Debye mass. The RGOPT remnant scale dependence at NLO is however significantly reduced in comparison: more precisely, for the largest shown (rescaled) coupling $g = 1$, the relative $P/P_0$ variation for $\pi T \le \mu \le 4\pi T$ is ≈8%, 1.5%, 0.8%, and 0.4%, respectively, for the $O(g^3)$ PT, SPT with screening PT mass, SPT with variational mass gap, and RGOPT.

## IV. RGOPT $\lambda\phi^4$ PRESSURE AT NNLO

Coming back to the basic free energy in Eq. (2.2), if neglecting the subtraction terms in Eq. (2.15), the formal lack of RG invariance from unmatched $m^4 \ln \mu$ terms remains relatively screened at one- and two-loop orders of thermal perturbative expansions for sufficiently small coupling, since perturbatively $m^4 \sim \lambda^2$. This can essentially explain why the remnant scale dependence of SPT remains quite moderate even at NLO, see Fig. 3, so that in comparison the NLO RGOPT improvement by merely a factor of 2 is not spectacular. But conversely, this can largely explain why a very sizeable scale dependence resurfaces for the NNLO SPT pressure [13], where the formally same order genuine three-loop $\lambda^2$ contributions are considered. In contrast, the RGOPT scale dependence is expected to further improve at higher orders, at least formally: being built on perturbatively restored RG invariance of the free energy at order $m^4\lambda^k$ for arbitrary $m$, the resulting mass gap exhibits a leading remnant scale dependence as $\tilde{m}^2 \sim \lambda T^2(1 + \cdots + O(\lambda^k \ln \mu))$. Thus, the dominant scale dependence in the free energy, coming from the leading term $\sim s_0 m^4 / \lambda$, is expected to appear first only at $O(\lambda^{k+1})$. Nevertheless, this expected trend could be largely spoiled, either by large perturbative coefficients (generically expected to grow at higher orders) or by the well-known thermal PT issues due to infrared divergent bosonic zero modes. It is thus important to investigate more explicitly the outcome of our construction at NNLO, where standard thermal PT starts to badly behave, to delineate the RGOPT scale dependence improvement that can be actually obtained. Moreover, concerning the $\lambda\phi^4$ model, the peculiar sign alternating beta-function coefficients $b_i$ from one- to three-loop orders, see Eq. (B3), implies that considering the running coupling alone, at the three-loop order, has a comparatively worse scale dependence than at the two-loop order, as illustrated in Fig. 4. This feature

<figure><img src="image_6.png" /><figcaption>FIG. 4. Relative scale dependence at successive orders of the running coupling in the $\lambda\phi^4$ model, given, respectively, at one-, two- and three-loop orders by Eqs. (3.8), (B8), and (B11), with $g \equiv \sqrt{\lambda/24}$ and $\pi T \le \mu \le 4\pi T$.</figcaption></figure>

tends to partly counteract the benefits of our RG-improved construction, when comparing NLO with NNLO.<fnref n="8" />

Applying the variational modification from Eqs. (3.1), (3.4) to the complete three-loop free energy, sum of Eq. (2.6) and Eq. (2.10), expanded consistently now to order $\delta^2$, and taking $\delta \to 1$, gives after algebra,

$$ (4\pi)^2 \mathcal{F}_0^{\delta^2} = \mathcal{E}_0^{\delta^2} - \frac{m^4}{8} \left( 3 + 2L + 8 \left( \frac{\gamma_0}{b_0} \right)^2 (L + J_2) \right) - \frac{T^4}{2} J_0 \\ + \frac{\gamma_0}{b_0} \left( \frac{\gamma_0}{b_0} - \frac{3}{2} \right) \frac{m^2 \Sigma_R}{\gamma_0 \lambda} + \frac{m^2}{32\pi^2 b_0} \Sigma_R (L + J_2) \\ + \frac{\Sigma_R^2}{128\pi^2 \gamma_0^2 \lambda} + F_{3l}, \qquad (4.1) $$

omitting the argument $m/T$ in thermal functions $J_i$, with $L \equiv \ln(\mu^2/m^2)$ and $\Sigma_R$ defined in Eq. (3.6). The original $O(\lambda^2)$ perturbative contribution, designated by $F_{3l}$ in Eq. (4.1), is already defined in Eq. (2.10), those terms being unaffected by Eq. (3.1) at the $\delta^2$ reexpansion order. The subtraction contributions in Eq. (4.1) after the modifications from Eq. (3.1) read explicitly, using $\gamma_0/b_0 = 1/6$,

$$ \mathcal{E}_0^{\delta^2} = -m^4 \left( \frac{14 s_0}{81 \lambda} + \frac{2}{9} s_1 + \frac{1}{3} s_2 \lambda + s_3 \lambda^2 \right). \quad (4.2) $$

The explicit expressions at NNLO for the RG and OPT Eqs. (3.3), (3.2) can be obtained straightforwardly from Eq. (4.1) after algebra. They are more involved than their NLO analogs in Eqs. (3.14), (3.15) and not particularly telling: some relevant expressions are given explicitly in Appendix A.

<footnote n="8">⁸We expect the corresponding behavior in QCD to be better, since the first three QCD beta-function coefficients have the same sign.</footnote>
Similarly to the NLO, at this stage without examining further constraints, one may *a priori* use any of the three possible (not independent) prescriptions to obtain the NNLO dressed optimized mass $\tilde{m}(\lambda, T)$: namely taking the solution of the OPT Eq. (3.2) or the massless RG Eq. (3.3) or the full RG Eq. (2.13). One may further exploit the freedom to incorporate the highest order subtraction term $s_3$ of Eq. (2.17) or not, the latter being formally a three-loop contribution but depending on four-loop RG coefficients, thus not necessary for NNLO RG invariance. The latter flexibility happens to give a relatively simple handle to circumvent the annoyance of possibly nonreal NNLO solutions that occur only at relatively large couplings in the $\lambda\phi^4$ model, quite similar to what happens at NLO. The behavior of those solutions for the different prescriptions is detailed in Appendix A for completeness.

As an outcome, in order to maximize the range of coupling and scale $\pi T \le \mu \le 4\pi T$ values where real solutions are obtained, it is appropriate to minimally neglect $s_3$ if using the OPT Eq. (3.2), and incorporating $s_3 \ne 0$ when using the RG Eq. (3.3). We mention that at NNLO, the full RG Eq. (2.13) gives no real solutions, at least for the relevant scale choice $\mu = 4\pi T$ that maximizes the values of $g(\mu)$ for a reference coupling $g \equiv g(2\pi T)$. [Actually one could recover real solutions only if truncating Eq. (2.13) maximally, keeping only $O(g^2)$ terms, but one then loses the crucial perturbative-matching properties, so that the corresponding solutions have to be rejected.] Some of these features could be intuitively expected: incorporating higher orders, either via $s_3 \ne 0$ or taking the more involved exact Eq. (2.13), renders the expression even more nonlinear in $m$, favoring the occurrence of nonreal solutions. But for the massless RG Eq. (3.3), provided that it has real solutions, incorporating $s_3$ could be expected to give better results, since the subtraction coefficients in Eq. (2.17) entering the pressure are originating directly from the RG coefficients.

Figure 5 illustrates as a function of the (rescaled) coupling $g$ our two different prescriptions thus retained at NNLO: respectively, applying to Eq. (4.1) the OPT Eq. (3.2) or the RG Eq. (3.3) [shown explicitly in Eq. (A1)], with resulting (unique) perturbatively matching exact solutions $\tilde{m}^{\text{OPT}}(\lambda, T)$ and $\tilde{m}^{\text{RG}}(\lambda, T)$, respectively. As is seen, despite being rather different functions of the coupling, they have similar very moderate scale dependence. Similarly to what happens at NLO, the RG solution is generically giving a slightly better scale dependence than the OPT one, since the former embeds more directly the perturbative RG properties. The previous exact $\tilde{m}$ solutions<fnref n="9" /> are then inserted within Eq. (4.1) to give the physical pressure. Figure 6 illustrates the corresponding pressures obtained from the two alternative OPT and RG mass

<figure><img src="image_7.png" /><figcaption>FIG. 5. NNLO (three-loops) $\tilde{m}$ from the OPT prescription Eq. (3.2) (with $s_3 = 0$) and RG prescription Eqs. (3.3) and (A1) (with $s_3 \ne 0$) as a function of $g \equiv \sqrt{\lambda/24}$, with the scale dependence $\pi T \le \mu \le 4\pi T$.</figcaption></figure>

prescriptions. As one can see, despite the quite different $\tilde{m}(g)$ in Fig. 5, the corresponding OPT and RG pressures are very close and similar in shape, and have comparable very moderate scale dependence. This feature appears as a very good indication of the previously mentioned expected convergence at a higher order of those *a priori* different prescriptions.

We illustrate our main results at successive LO, NLO, and NNLO in Fig. 7, compared with the PT pressure, Eq. (2.12), and SPT pressure at successive orders. The NNLO PT, with successive terms up to $O(\lambda^2)$, has a substantially larger remnant scale dependence than the $O(\lambda^{3/2})$ PT in Fig. 3. Concerning the SPT, at NNLO, we only show the results from using the mass gap Eq. (3.19), following the very same prescription as in [13], namely using also Eq. (3.19) at NNLO. [N.B. using the perturbative screening mass instead, Eq. (3.20), gives a much larger scale dependence, that we do not illustrate.] In contrast, one can see that the RGOPT pressure is remarkably stable from comparing LO to NLO and NNLO, and has a very

<figure><img src="image_8.png" /><figcaption>FIG. 6. Comparison of NNLO pressure $P/P_{\text{ideal}}$ for the OPT and RG $\tilde{m}$ prescriptions (with $s_3 = 0$ and $s_3 \ne 0$, respectively), as a function of $g \equiv \sqrt{\lambda/24}$ with $\pi T \le \mu \le 4\pi T$. NB: the lower pressure values correspond to the largest $\mu = 4\pi T$ values due to infrared freedom of $\lambda\phi^4$.</figcaption></figure>

<footnote n="9">⁹At NNLO, we use the exact expressions of all thermal integrals, in particular, Eqs. (C1), (C4), as their high-$T$ expressions in Eqs. (C7), (C8) are not very good approximations.</footnote>
<figure><img src="image_9.png" /><figcaption>FIG. 7. RGOPT pressure at successive LO,NLO,NNLO orders versus NLO,NNLO PT and NLO,NNLO SPT pressures, with ($g \equiv \sqrt{\lambda}/24$), $\pi T \le \mu \le 4\pi T$. NB: the lower pressure values correspond to the largest $\mu = 4\pi T$ values due to infrared freedom of $\lambda\phi^4$.</figcaption></figure>

moderate remnant scale dependence at NNLO, almost invisible at the figure scale until a relatively large $g \approx 0.8$. Actually, the improvement from NLO to NNLO becomes only moderate for relatively large coupling values $0.9 \lesssim g \lesssim 1$, that we understand as the countereffect from the worse scale dependence of the sole NNLO running coupling as above explained; see Fig. 4. Overall, the scale dependence is drastically improved as compared to PT and SPT.

We remark finally that in principle a yet more involved prescription could be to combine Eqs. (3.2) and (3.3), thus determining variational dressed mass and coupling, as was indeed investigated at zero temperature in [25] and briefly explored at finite temperature at NLO in [28]. However, due to the much more involved (very nonlinear) $m, T$ dependences, besides the numerical complications to simultaneously solve those nonlinear equations at finite $T$, one issue is that the range where both equations have a (common) real solution happens to be much more restricted. To extend this range would need at least to exploit more complicated prescriptions with renormalization scheme changes introducing extra variational parameters, generalizing what was done for QCD in [31,33], that is beyond our present scope.

## V. CONCLUSIONS AND OUTLOOK

We have illustrated in the $\lambda\phi^4$ model up to NNLO how the RGOPT resummation of thermal perturbative expansions, consistently maintaining perturbative RG invariance, leads to drastically improved convergence and remnant renormalization scale dependence at successive orders, as

compared with PT and with related thermal perturbative resummation approaches such as SPT. In particular, RGOPT remains very stable from NLO to NNLO up to relatively large coupling values. We have compared two different prescriptions a priori available in our framework to defining a resummed dressed mass at a given order, that give very close and stable results. As could be intuitively expected, the prescription defining the dressed mass from the RG Eq. (3.3), thus more directly embedding RG properties, appears to give a slightly better remnant scale dependence than the more traditional variational mass prescription from Eq. (3.2) optimizing the pressure.

The RGOPT has been applied recently at NLO in thermal QCD [32,33], yet only in the quark sector, thus treating the gluon contributions apart purely perturbatively, due to some present technical limitations. Although such a relatively simple approximation shows a very good agreement with lattice simulations down to relatively low temperatures near the pseudotransition, it is clearly an important next step to extend our approach to similarly treat the crucial gluon sector, largely responsible for the poorly convergent weak coupling expansion of thermal QCD, due to zero mode infrared divergences that start to show up at NNLO. Since the very stable NNLO RGOPT properties here demonstrated for the scalar model entails an appropriate resummation to all orders of infrared divergent bosonic zero modes, we anticipate similarly good properties to hold also in the QCD gluon sector, once the technical (computational) difficulties to readily adapt the RGOPT approach to the gluon sector will be overcome. More precisely, the RG-restoring subtraction contributions analogous to Eq. (2.15), necessarily involve gluon mass terms. The well-known explicitly gauge-invariant framework that entails, among other features, a screening gluon mass term consistent at finite temperature, is the HTL nonlocal Lagrangian formalism [20]. However, the contributions analogous to Eq. (2.15) require at NLO some very involved HTL integral calculations, not yet fully available in the literature, that we leave for future work.

## APPENDIX A: PROPERTIES OF RG AND OPT SOLUTIONS AT NNLO

This Appendix examines in some details the properties of the different possible prescriptions at NNLO. We first give the explicit expression of the NNLO RG Eq. (3.3), straightforward to obtain after algebra from Eq. (4.1) that takes the rather compact form,

$$f_{\text{RG}}^{\text{NNLO}} = -\frac{m^4}{2} \left( 1 - 6 \frac{\gamma_0}{b_0} + 8 \left( \frac{\gamma_0}{b_0} \right)^2 \right) + \frac{m^2}{16\pi^2 b_0} \left( \Sigma_R - m^2 \gamma_0 \lambda (L + J_2) \right) - \frac{m^2 \Sigma_R}{32\pi^2 \gamma_0} + G_{3l} \\ + \beta^{(3)}(\lambda) \left\{ \frac{m^2 \Sigma_R}{32\pi^2 b_0 \lambda} (L + J_2) + \frac{\Sigma_R^2}{128\pi^2 \gamma_0^2 \lambda^2} + \frac{2}{\lambda} F_{3l} - m^4 \left( -\frac{14 s_0}{81 \lambda^2} + \frac{s_2}{3} + 2s_3 \lambda \right) \right\} = 0, \quad (\text{A1})$$
where

$$G_{3l} = -\frac{1}{24} \left(\frac{\lambda}{16\pi^2}\right)^2 \left[m^4 \left(15L^2 + 34L + \frac{41}{2} + 6(L+1)J_2\right) - m^2 T^2 J_1(24L + 28 + 6J_2) + 9T^4 J_1^2\right], \quad (A2)$$

$$\beta^{(3)}(\lambda) = \lambda^2(b_0 + b_1\lambda + b_2\lambda^2), \qquad (A3)$$

with $F_{3l}$ given in Eq. (2.10) and $\Sigma_R$ in Eq. (3.6). The OPT Eq. (3.2) at NNLO is similarly easily obtained from Eq. (4.1), but it gives a somewhat more lengthy expression that we thus refrain to give explicitly.

Figure 8 illustrates, for $s_3 = 0$ (top) and $s_3 \neq 0$ (bottom), the behavior of these exact RG Eq. (A1) and OPT (3.2), respectively, as function of $m$ for two representative $g \equiv \sqrt{\lambda/24}$ values and for the renormalization scale $\mu = 4\pi T$, giving the largest coupling for a given $g$ input, thus the most problematic case for sufficiently large $g$. For the simpler $s_3 = 0$ prescription, the RG equation has a unique perturbative-matching solution for $0 < g \lesssim .75$ (for $\mu = 4\pi T$),

<figure><img src="image_10.png" /><figcaption>FIG. 8. Top panel: OPT (dashed lines) and RG (thick lines) equations at NNLO as a function of $m$, for $s_3 = 0$, taking three-loop order running coupling $g(4\pi T)$ ($g \equiv \sqrt{\lambda/24}$) for two different values of the reference coupling, $g(2\pi T) = 0.6$ and $g(2\pi T) = 0.8$, respectively. Bottom panel: same captions for $s_3 \neq 0$.</figcaption></figure>

but above this maximal $g$ value, an inflection point appears that pushes the previous perturbative-matching solution to complex values, although being close to real values as can be seen in Fig. 8 (top). Note also that another RG solution appears at a higher mass value, but the latter is not matching standard perturbative behavior, so it has to be rejected. Rather similar features are obtained for $\mu = 2\pi T$, but the disappearance of the real perturbative-matching solution is delayed to higher $g \gtrsim 0.95$ values. In contrast for $s_3 \neq 0$, a similar behavior is obtained but the disappearance of a real perturbative-matching RG solution is delayed to substantially higher $g \gtrsim 1$ values for any $\mu \le 4\pi T$: accordingly, for $0 < g \le 1$ and $\pi T \le \mu \le 4\pi T$ the perturbative-matching solution is real and unique. The OPT equation has a unique perturbative-matching solution for both $s_3 = 0$ or $s_3 \neq 0$, and we adopt the simpler minimal choice $s_3 = 0$ for our NNLO results. Notice that despite their very different expressions, the RG and OPT equations give $\tilde{m}$ solutions that are rather close to each other. (Even when the RG solution is not real, it is seen to be very close to the OPT real one.)

Finally, note that the full RG Eq. (2.13) has no real solutions for $0 < g \le 1$ for $\mu = 4\pi T$: the latter is not illustrated on Fig. 8, but its shape looks quite similar to the RG one in the range where the latter gives complex solutions (also being close to give a real solution).

## APPENDIX B: RG INGREDIENTS AND COUNTERTERMS

Our normalization for the beta function and mass anomalous dimensions are, respectively,

$$\beta(\lambda) \equiv \frac{d\lambda}{d\ln\mu} = b_0\lambda^2 + b_1\lambda^3 + b_2\lambda^4 + \dots \quad (B1)$$

and

$$\gamma_m(\lambda) \equiv \frac{d \ln m}{d \ln \mu} = \gamma_0 \lambda + \gamma_1 \lambda^2 + \gamma_2 \lambda^3 + \dots, \quad (B2)$$

where up to the three-loop order [39],
$$ (4\pi)^2 b_0 = 3; \quad (4\pi)^4 b_1 = -17/3; \quad (4\pi)^6 b_2 = \frac{3915 + 2592\zeta(3)}{216}, \qquad (B3) $$

$$ (4\pi)^2 \gamma_0 = 1/2; \quad (4\pi)^4 \gamma_1 = -5/12; \quad (4\pi)^6 \gamma_2 = \frac{7}{4}. \qquad (B4) $$

Next, the coupling and mass counterterm are expressed as perturbative $\lambda$ and $1/\epsilon$ expansions, with ($d=4-2\epsilon$), $\lambda_0 \equiv \mu^{2\epsilon}\lambda Z_\lambda$, $m_0 \equiv mZ_m$ (where only the $\lambda^2$ terms are needed at the three-loop order),

$$ Z_{\lambda} = 1 + \frac{b_0}{2\epsilon}\lambda + \left[ \left(\frac{b_0}{2\epsilon}\right)^2 + \frac{b_1}{4\epsilon} \right] \lambda^2 + \mathcal{O}(\lambda^3) \qquad (B5) $$

$$ Z_m = 1 + \frac{\gamma_0}{2\epsilon}\lambda + \left[ \frac{\gamma_0(\gamma_0 + b_0)}{8\epsilon^2} + \frac{\gamma_1}{4\epsilon} \right] \lambda^2 + \mathcal{O}(\lambda^3). \qquad (B6) $$

The necessary additional vacuum energy counterterms up to three loop order read [37] in our normalization conventions,

$$ (4\pi)^2 \Delta\mathcal{F}_0 = m^4 \left[ \frac{1}{4\epsilon} + \left(\frac{\lambda}{16\pi^2}\right) \frac{1}{8\epsilon^2} + \left(\frac{\lambda}{16\pi^2}\right)^2 \left( \frac{5}{48\epsilon^3} - \frac{5}{72\epsilon^2} + \frac{1}{96\epsilon} \right) \right]. \qquad (B7) $$

Next, the (exact) two-loop running coupling used in the numerics is

$$ \lambda^{(2-\text{loop})}(\mu) = \frac{\lambda(\mu_0)}{f_W(\lambda(\mu_0), \text{ln}\frac{\mu}{\mu_0})}, \qquad (B8) $$

with

$$ f_W(\lambda, L_\mu) = 1 - b_0 L_\mu \lambda + \frac{b_1}{b_0} \lambda \ln \left( f_W \frac{1 + \frac{b_1}{b_0} \lambda f_W^{-1}}{1 + \frac{b_1}{b_0} \lambda} \right) = - \frac{b_1}{b_0} \lambda \left\{ 1 + W \left[ - \left( 1 + \frac{b_0}{b_1 \lambda} \right) e^{-\left[ 1 + \frac{b_0}{b_1 \lambda} (1 - b_0 \lambda L) \right]} \right] \right\}, \qquad (B9) $$

where $L_\mu \equiv \ln(\mu/\mu_0)$ and $W(x) \equiv \ln(W/x)$ is the Lambert implicit function. For the range of coupling values illustrated in our main figures, $g \equiv \sqrt{\lambda/24} \lesssim 1$, Eqs. (B8), (B9) do not give much visible differences with a simpler perturbatively truncated expansion at order $\lambda^3$,

$$ \lambda^{-1}(\mu) \approx \lambda^{-1}(\mu_0) - b_0 L_\mu - (b_1 L_\mu)\lambda - \left(\frac{1}{2}b_0 b_1 L_\mu^2\right)\lambda^2 - \left(\frac{1}{2}b_1^2 L_\mu^2 + \frac{1}{3}b_0^2 b_1 L_\mu^3\right)\lambda^3 + \mathcal{O}(\lambda^4). \qquad (B10) $$

At the three-loop order, we used quite similarly an exact integral giving the running coupling $\lambda(\mu)$ as a more involved implicit function, numerically solved for $\mu$ as a function of $\mu_0$. For coupling values not too large, there is not much visible differences with a more common perturbatively truncated running coupling,

$$ \lambda^{(3-\text{loop})}(\mu) = \frac{\lambda(\mu_0)}{f^{3L}}, \\ f^{3L} = 1 - b_0 L_\mu \lambda - b_1 L_\mu \lambda^2 - \left(b_0 \frac{b_1}{2} L_\mu^2 + b_2 L_\mu\right) \lambda^3 \\ + \frac{L_\mu}{6b_0^2} \left(6b_1(b_1^2 - 2b_0b_2) - 3b_0^2(b_1^2 + 2b_0b_2)L_\mu - 2b_0^4 b_1 L_\mu^2\right) \lambda^4. \qquad (B11) $$

## APPENDIX C: TWO- AND THREE-LOOP IRREDUCIBLE INTEGRALS

For completeness, we reproduce here the explicit expressions of the thermal three-loop massive integrals $K_2$ and $K_3$, originally calculated in [34], and entering Eq. (2.10),
$$ K_2\left(\frac{m}{T}\right) = -\frac{32}{T^4} \int_0^\infty dp \, p \frac{n(E_p)}{E_p} \int_0^p dq \, q \frac{n(E_q)}{E_q} \int_{p-q}^{p+q} dk \, k \sum_{\sigma=-1,+1} f_2(E_\sigma, k), \quad (C1) $$

where $E_p = (p^2 + m^2)^{1/2}$, $n(x) = (e^{x/T} - 1)^{-1}$,

$$ \begin{aligned} f_2(E, k) &= \left(\frac{E^2 - M_k^2}{E^2 - k^2}\right)^{\frac{1}{2}} \ln \frac{(E^2 - k^2)^{\frac{1}{2}} + (E^2 - M_k^2)^{\frac{1}{2}}}{(E^2 - k^2)^{\frac{1}{2}} - (E^2 - M_k^2)^{\frac{1}{2}}}, && k^2 < E^2 - 4m^2 \\ &= 2 \left(\frac{M_k^2 - E^2}{E^2 - k^2}\right)^{\frac{1}{2}} \arctan \left(\frac{E^2 - k^2}{M_k^2 - E^2}\right)^{\frac{1}{2}}, && E^2 - 4m^2 < k^2 < E^2 \\ &= \left(\frac{M_k^2 - E^2}{k^2 - E^2}\right)^{\frac{1}{2}} \ln \frac{(M_k^2 - E^2)^{\frac{1}{2}} + (k^2 - E^2)^{\frac{1}{2}}}{(M_k^2 - E^2)^{\frac{1}{2}} - (k^2 - E^2)^{\frac{1}{2}}}, && E^2 < k^2 \end{aligned} \quad (C2) $$

and

$$ \begin{aligned} M_k^2 &= 4m^2 + k^2, && k \equiv |\mathbf{p} + \mathbf{q}|, \\ E_\sigma(p, q) &= \sqrt{p^2 + m^2 + \sigma\sqrt{q^2 + m^2}}. \end{aligned} \quad (C3) $$

$$ K_3\left(\frac{m}{T}\right) = \frac{96}{T^4} \int_0^\infty dp \, p \frac{n(E_p)}{E_p} \int_0^p dq \, q \frac{n(E_q)}{E_q} \int_0^q dr \, r \frac{n(E_r)}{E_r} \\ \times \sum_{\sigma, \tau = -1, +1} \{f_3(E_{\sigma\tau}, p+q+r) - f_3(E_{\sigma\tau}, p+q-r) - f_3(E_{\sigma\tau}, p-q+r) + f_3(E_{\sigma\tau}, p-q-r)\}, \quad (C4) $$

where

$$ \begin{aligned} f_3(E, p) &= p \ln \frac{m^2 - E^2 + p^2}{m^2} + 2(m^2 - E^2)^{\frac{1}{2}} \arctan \frac{p}{(m^2 - E^2)^{\frac{1}{2}}}, && E^2 < m^2 \\ &= p \ln \frac{|E^2 - m^2 - p^2|}{m^2} + (E^2 - m^2)^{\frac{1}{2}} \ln \frac{(E^2 - m^2)^{\frac{1}{2}} + p}{|(E^2 - m^2)^{\frac{1}{2}} - p|}, && E^2 > m^2 \end{aligned} \quad (C5) $$

and

$$ E_{\sigma\tau}(p, q, r) = \sqrt{p^2 + m^2} + \sigma\sqrt{q^2 + m^2} + \tau\sqrt{r^2 + m^2}. \quad (C6) $$

In the limit $x \equiv m/T \to 0$, Eqs. (C1), (C4) can be expressed analytically as

$$ K_2(x) \approx \frac{(4\pi)^4}{72} \left( \ln x + \frac{1}{2} + \frac{\zeta'(-1)}{\zeta(-1)} \right) - 372.65x(\ln x + 1.4658), \quad (C7) $$

$$ K_3(x) \approx \frac{(4\pi)^4}{48} \left( -\frac{7}{15} + \frac{\zeta'(-1)}{\zeta(-1)} - \frac{\zeta'(-3)}{\zeta(-3)} \right) + 1600.0x(\ln x + 1.3045). \quad (C8) $$

For the numerical results, we rather use the exact expressions Eqs. (C1), (C4).
[1] J. P. Blaizot, E. Iancu, and A. Rebhan, *Thermodynamics of the High Temperature Quark Gluon Plasma*, edited by R. C. Hwa et al., Quark gluon plasma (World Scientific, Singapore, 2004), pp. 60–122; U. Kraemmer and A. Rebhan, *Rep. Prog. Phys.* **67**, 351 (2004).

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

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

[4] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz, and K. K. Szabo, *Nature (London)* **443**, 675 (2006); Y. Aoki, Z. Fodor, S. D. Katz, and K. K. Szabo, *Phys. Lett. B* **643**, 46 (2006); M. Cheng et al., *Phys. Rev. D* **74**, 054507 (2006).

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

[6] R. R. Parwani, *Phys. Rev. D* **45**, 4695 (1992); **48**, 5965(E) (1993).

[7] F. Karsch, A. Patkos, and P. Petreczky, *Phys. Lett. B* **401**, 69 (1997).

[8] J. P. Blaizot, A. Ipp, and N. Wschebor, *Nucl. Phys. A* **849**, 165 (2011).

[9] J. P. Blaizot, E. Iancu, and A. Rebhan, *Phys. Rev. D* **63**, 065003 (2001).

[10] J. Berges, S. Borsanyi, U. Reinosa, and J. Serreau, *Phys. Rev. D* **71**, 105004 (2005); Ann. Phys. (Amsterdam) **320**, 344 (2005).

[11] M. E. Carrington, B. A. Meggison, and D. Pickering, *Phys. Rev. D* **94**, 025018 (2016).

[12] J. P. Blaizot and N. Wschebor, *Phys. Lett. B* **741**, 310 (2015).

[13] J. O. Andersen, E. Braaten, and M. Strickland, *Phys. Rev. D* **63**, 105008 (2001).

[14] J. O. Andersen and M. Strickland, *Phys. Rev. D* **64**, 105012 (2001); Ann. Phys. (Amsterdam) **317**, 281 (2005).

[15] R. L. S. Farias, G. Krein, and R. O. Ramos, *Phys. Rev. D* **78**, 065046 (2008).

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

[17] V. I. Yukalov, *Theor. Math. Phys.* **28**, 652 (1976); W. E. Caswell, *Ann. Phys. (N.Y.)* **123**, 153 (1979); I. G. Halliday and P. Suranyi, *Phys. Lett. B* **85B**, 421 (1979); R. P. Feynman and H. Kleinert, *Phys. Rev. A* **34**, 5080 (1986); A. Duncan and M. Moshe, *Phys. Lett. B* **215**, 352 (1988); H. F. Jones and M. Moshe, *Phys. Lett. B* **234**, 492 (1990); A. Neveu, *Nucl. Phys. B*, *Proc. Suppl.* **18**, 242 (1991); V. Yukalov, *J. Math. Phys (N.Y.)* **32**, 1235 (1991); S. Gandhi, H. F. Jones, and M. Pinto, *Nucl. Phys. B* **359**, 429 (1991); C. M. Bender, F. Cooper, K. A. Milton, M. Moshe, S. S. Pinsky, and L. M. Simmons, *Phys. Rev. D* **45**, 1248 (1992); H. Yamada, *Z. Phys. C* **59**, 67 (1993); A. N. Sisakian, I. L. Solovtsov, and

O. P. Solovtsova, *Phys. Lett. B* **321**, 381 (1994); H. Kleinert, *Phys. Rev. D* **57**, 2264 (1998); *Phys. Lett. B* **434**, 74 (1998).

[18] R. Seznec and J. Zinn-Justin, *J. Math. Phys. (N.Y.)* **20**, 1398 (1979); J. C. Le Guillou and J. Zinn-Justin, *Ann. Phys. (N.Y.)* **147**, 57 (1983).

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

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

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

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

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

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

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

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

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

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

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

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

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

[32] J. L. Kneur, M. B. Pinto, and T. E. Restrepo, *Phys. Rev. D* **104**, L031502 (2021).

[33] J. L. Kneur, M. B. Pinto, and T. E. Restrepo, *Phys. Rev. D* **104**, 034003 (2021).

[34] J. O. Andersen, E. Braaten, and M. Strickland, *Phys. Rev. D* **62**, 045004 (2000).

[35] J. Frenkel, A. V. Saa, and J. C. Taylor, *Phys. Rev. D* **46**, 3670 (1992); P. B. Arnold and C. X. Zhai, *Phys. Rev. D* **50**, 7603 (1994); R. Parwani and H. Singh, *Phys. Rev. D* **51**, 4518 (1995).

[36] E. Braaten and A. Nieto, *Phys. Rev. D* **51**, 6990 (1995).

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

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

[39] E. Brezin, J. C. Le Guillou, and J. Zinn-Justin, *Phys. Rev. D* **9**, 1121 (1974).