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

Loïc Fernandez  Laboratoire Charles Coulomb (L2C), UMR 5221 CNRS-Univ.
Montpellier, Montpellier, France  Jean-Loïc Kneur  Laboratoire Charles Coulomb
(L2C), UMR 5221 CNRS-Univ. Montpellier, Montpellier, France

###### Abstract

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.

##  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. Trev ; kapusta-gale ; la-vu  for reviews), it leads to nonanalytical
terms in the coupling, that happen to give poorly convergent successive
orders, even when pushed to the highest perturbative order available,
furthermore with increasingly sizeable remnant scale dependence. This
situation is notoriously illustrated for asymptotically free thermal QCD,
where however lattice simulations (LS) offer a powerful genuinely
nonperturbative alternative 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 QCDlattT . Nevertheless, the famous numerical sign problem
signpb  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\geq 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)SPTearly ; SPT , the
nonperturbative renormalization group (NPRG)NPRG  approach, the two-particle
irreducible (2PI) formalism2PI ; 2PI3loop ; 2PI4loop , or other
approachesBlaizot-Wschebor . 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-loopSPT3l ; SPT3lnext ; OPT3l  and even four-loop
ordersspt4L . SPT may also be viewed as a particular case, in the thermal
context, of the so-called optimized perturbation theory (OPT)111OPT and its
many variants appear under different names in the literature OPT-LDE ; odm ;
pms ., 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 resummation of perturbative
expansion. The generalization of SPT tailored to treat the much more involved
thermal gauge theoriesHTL , the Hard Thermal Loop perturbation theory
(HTLpt)htlpt1 , has been pushed to three-loop order HTLPT3loop ; HTLPTMU . The
three-loop resultsHTLPTMU  agreement with LS is quite remarkable down to about
twice the critical temperature, for a renormalization scale choice $\sim 2\pi
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)rgopt1 ;
rgopt_alphas . 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 rgopt_alphas  of the basic
QCD scale ($\Lambda_{\overline{\rm MS}}$) or related coupling $\alpha_{S}$,
and of the quark condensate rgopt_qq ; rgopt_qq5 . The RGOPT was extended at
finite temperatures for the scalar $\lambda\phi^{4}$ in prl_phi4 ; prd_phi4
and for the non-linear sigma model (NLSM) in nlsm , 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 temperaturesprdCOLD ,
and at finite temperature and densityrgopt_NLOQCD1 ; rgopt_NLOQCD2 , 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 $\lambda\phi^{4}$ 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 $\lambda\phi^{4}$ 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 prl_phi4 ; prd_phi4 . 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.

##  II Review of standard (massive) thermal perturbative expansion

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

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

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

Figure 1:  Free energy diagrams up to NNLO in $\lambda\phi^{4}$ model.

###  II.1 Free energy up to NNLO

We first recall the resultkapusta-gale ; SPT3l  for the two-loop free energy
(equivalently minus the pressure) including a mass term, corresponding to the
first two graphs 222In Fig.1 counterterm graphs are omitted for simplicity. in
Fig.1,

${\cal F}_{0}=\frac{1}{2}\hbox{$\sum$}\\!\\!\\!\\!\\!\\!\\!\int_{P}\ln({\bf
P}^{2}+m_{0}^{2})+\frac{\lambda_{0}}{8}\left(\hbox{$\sum$}\\!\\!\\!\\!\\!\\!\\!\int_{P}\frac{1}{{\bf
P}^{2}+m_{0}^{2}}\right)^{2}+{\cal F}_{0}^{2l,\rm ct},$ (2)

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

$\hbox{$\sum$}\\!\\!\\!\\!\\!\\!\\!\int_{P}\equiv
T\sum_{n}(\frac{\mu^{2}e^{\gamma_{E}}}{4\pi})^{\epsilon}\,\int\frac{d^{3-2\epsilon}{\bf
p}}{(2\pi)^{3-2\epsilon}}.$ (3)

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

$\displaystyle{\cal
F}^{3l}_{0}=-\frac{\lambda_{0}^{2}}{48}\left[3\left(\hbox{$\sum$}\\!\\!\\!\\!\\!\\!\\!\int_{p}\frac{1}{{\bf
p}^{2}+m_{0}^{2}}\right)^{2}\hbox{$\sum$}\\!\\!\\!\\!\\!\\!\\!\int_{q}\frac{1}{({\bf
q}^{2}+m_{0}^{2})^{2}}\right.$
$\displaystyle\left.+\hbox{$\sum$}\\!\\!\\!\\!\\!\\!\\!\int_{pqr}\frac{1}{({\bf
p}^{2}+m_{0}^{2})({\bf q}^{2}+m_{0}^{2})({\bf r}^{2}+m_{0}^{2})({\bf
p+q+r}^{2}+m_{0}^{2})}\right]$ $\displaystyle+{\cal F}_{0}^{3l,\rm ct}.$ (4)

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

$\lambda_{0}=\lambda Z_{\lambda},\;\;m_{0}=m\,Z_{m},$ (5)

where $Z_{\lambda},Z_{m}$ are the standard coupling and mass counterterms for
the massive $\lambda\phi^{4}$ model, given for completeness respectively in
Eqs.(47),(48) in Appendix B. Those are expanded for the relevant three-loop
order case here to perturbative order $\lambda^{2}$. The remaining divergences
are then removed by an additive renormalization from the vacuum energy
countertermsSPT3l , formally represented as ${\cal F}_{0}^{2l,\rm ct}$, ${\cal
F}_{0}^{3l,\rm ct}$ in the above expressions, and also given for completeness
in Eq.(49) in Appendix B.

Once the mass, coupling, and vacuum energy counterterms have been accounted to
cancel the original divergences, one obtains the ($\overline{\rm MS}$-scheme)
renormalized free energy. The one- and two-loop contributions read SPT3l ;
kapusta-gale :

$(4\pi)^{2}{\cal F}^{NLO}={\cal
E}_{0}-\frac{1}{8}m^{4}\left(3+2L\right)-\frac{T^{4}}{2}J_{0}(\frac{m}{T})+\frac{1}{8}(\frac{\lambda}{16\pi^{2}})\left[(L+1)m^{2}-T^{2}J_{1}(\frac{m}{T})\right]^{2}\,,$
(6)

where $L\equiv\ln(\mu^{2}/m^{2})$ and we explicitly separated the thermal and
non-thermal contributions. Here and in all related renormalized expressions
below, $\mu$ stands for the arbitrary renormalization scale introduced by
dimensional regularization in the $\overline{\rm MS}$-scheme, Eq.(3), and
$\lambda\equiv\lambda(\mu)$. Note carefully that ${\cal E}_{0}$ in Eq.(6)
represents a finite vacuum energy term, to be speficied below, that plays a
crucial role in our approach as will be reexamined in Sec. II.2.
The standard (dimensionless) thermal integrals appearing in Eq.(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}}}\>\frac{1}{e^{\sqrt{t^{2}+x^{2}}}-1}\,,$
(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)\,.$ (8)

Also, a high-$T$ expansionkapusta-gale  such as

$J_{0}(x)\simeq\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]+{\cal
O}(x^{6})\,,$ (9)

is often a rather good approximation as long as
$x\raisebox{-3.69899pt}{~{}\shortstack{$<$ \\\ [-0.07cm] $\sim$}}~{}1$, i.e.,
$m\ll T$.

Next the NNLO contribution involves the additional genuine massive three-loop
second integral in Eq.(4), first calculated in scal3l . After algebra, the
complete three-loop contribution can be expressed asSPT3l

$\displaystyle
F_{3l}=-\frac{1}{48}(\frac{\lambda}{16\pi^{2}})^{2}\left[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}(\frac{m}{T})\right)\right.$
$\displaystyle\left.-m^{2}T^{2}J_{1}(\frac{m}{T})\left(12L^{2}+28L-12-\pi^{2}-4c_{1}+6(L+1)J_{2}(\frac{m}{T})\right)\right.$
$\displaystyle\left.+T^{4}\left(3(3L+4+J_{2}(\frac{m}{T}))J^{2}_{1}(\frac{m}{T})+6K_{2}(\frac{m}{T})+4K_{3}(\frac{m}{T})\right)\right]$
(10)

where

$c_{0}=\frac{275}{12}+\frac{23}{2}\zeta(2)-2\zeta(3)\simeq
39.429,\;\;c_{1}=-\frac{59}{8}-\frac{3}{2}\zeta(2)\simeq-9.8424,$ (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.scal3l , reproduced for
selfcontainedness in Eqs.(54),(57) in Appendix C.

To complete this subsection on basic thermal perturbative expressions, for
latter easier reference and comparison purpose we also recall the expression
of the (massless) PT pressurePptphi4 ; mdeb ; Trev  up to NNLO:

$\frac{P}{P_{0}}=1-\frac{5}{4}(\frac{\lambda}{16\pi^{2}})+5\frac{\sqrt{6}}{3}(\frac{\lambda}{16\pi^{2}})^{3/2}+\frac{15}{4}(\frac{\lambda}{16\pi^{2}})^{2}\,(\ln\frac{\mu}{2\pi
T}+0.40)+{\cal O}(\lambda^{5/2}),$ (12)

where $P_{0}=(\pi^{2}/90)T^{4}$ is the ideal bosonic gas pressure.

###  II.2 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.(6), with ${\cal 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}\;,$ (13)

with $\beta(\lambda)$, $\gamma_{m}(\lambda)$ given in Eqs.(43),(44), 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 lowest orders
terms from $\beta(\lambda)$ or $\gamma_{m}(\lambda)$ in Eq.(13), those being
at least of next order ${\cal O}(\lambda m^{4})$. This is a manifestation that
perturbative RG invariance generally occurs from cancellations between terms
from the RG equation at 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 ${\cal E}_{0}$ to
the action, although this term is usually ignored, minimally set to zero in
the (thermal) literature SPT3l ; htlpt1 ; HTLPT3loop . Following Refs.
prl_phi4 ; prd_phi4  the easiest way to construct a (perturbatively) RG-
invariant renormalized vacuum energy is to determine ${\cal E}_{0}$ order by
order as a perturbative series from the reminder of Eq. (13) applied on the
non RG-invariant finite part of Eq. (6):

$\mu\frac{d}{d\mu}{\rm{\cal E}_{0}}(\lambda,m)\equiv-{\rm
Remnant}(\lambda,m)=-\mu\frac{d}{d\mu}[{\cal F}_{0}({\rm{\cal E}_{0}\equiv
0})|_{\rm finite}]\;.$ (14)

Accordingly ${\cal E}_{0}$ has the form

${\rm{\cal E}_{0}}(\lambda,m)=-\frac{m^{4}}{\lambda}\sum_{k\geq
0}s_{k}\lambda^{k}\;,$ (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) rgopt_alphas .
This procedure leaves non RG-invariant remnant terms of higher orders, that
may be treated similarly once higher order terms are considered. Explicitly,
we obtainprd_phi4

$\displaystyle s_{0}=\frac{1}{2(b_{0}-4\gamma_{0})}=8\pi^{2}\;$ $\displaystyle
s_{1}=\frac{(b_{1}-4\gamma_{1})}{8\gamma_{0}\,(b_{0}-4\gamma_{0})}=-1\;$ (16)

and similarly for the next orders presently relevant,

$\displaystyle
s_{2}=\frac{96\pi^{2}\,(b_{0}-128\pi^{2}\left((1+4s_{1})\gamma_{1}-s_{0}(b_{2}-4\gamma_{2})\right)-41}{12288\,\pi^{4}(b_{0}+4\gamma_{0})}=\frac{23+36\zeta[3]}{480\,\pi^{2}}\simeq
0.01399\;,$ $\displaystyle
s_{3}=\frac{-709+12\pi^{4}-2628\zeta(3)-5400\zeta(5)}{720\,(16\pi^{2})^{2}}.$
(17)

The explicit RG coefficients in intermediate expressions emphasizes 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.(14) suggests, its
RHS precisely defines the vacuum energy anomalous dimension, that has been
calculated even to five-loop order for the general $O(N)$ scalar model
vacanom_kastening . Our independent results for the $s_{k}$ are fully
consistent with vacanom_kastening . A subtlety is that according to Eq. (15),
$s_{k}$ is strictly required for perturbative RG invariance at order
$\lambda^{k}$, but contributes at order $\lambda^{k-1}$. So at order
$\lambda^{k}$ one may choose minimally to include only $s_{0},\cdots s_{k}$,
or more completely include also $s_{k+1}$, incorporating in this way 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 prl_phi4 ; prd_phi4  at LO and NLO, before to
extend our approach at the technically more involved NNLO. We also underline
some important features that were not plainly discussed in prd_phi4 . After
restoring in a first stage perturbative RG invariance of the massive free
energy, leading to the crucial additional term in Eq.(15), one performs on the
resulting complete expression the variational modification, according to

${\cal F}_{0}(m^{2}\to m^{2}\,(1-\delta)^{2a}\,,\;\;\lambda\to\delta\lambda),$
(18)

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.(18) at
successive orders, $\delta^{k}$, at the same order as the original
perturbative expression, and set $\delta\to 1$ afterwards. This leaves a
remnant $m$ dependence at any order $k$, that may be conveniently fixed by a
stationarity prescription pms ,

$\frac{\partial{\cal F}_{0}^{(k)}}{\partial
m}(m,\lambda,\delta=1)|_{m\equiv\tilde{m}}\equiv 0,$ (19)

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

$\left[\mu\frac{\partial}{\partial\mu}+\beta(\lambda)\frac{\partial}{\partial\lambda}\right]{\cal
F}_{0}^{(k)}(m,\lambda,a,\delta=1)=0,$ (20)

that at leading RG order uniquely fixesrgopt_alphas ; prd_phi4

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

simply in terms of the universal (renormalization scheme independent) first
order RG coefficients333At 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.(21). But this
would involve extra arbitrary variational parameters with no compelling
reasons. We thus keep the simpler form Eq.(21) at successive orders for
sensible comparisons..

At higher orders, Eq.(20) 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.(19). Importantly Eq.(21) guarantees in
addition that either Eq.(19) or Eq.(20) have at least one (essentially unique)
solution matching the $T=0$ perturbative behavior rgopt_alphas ; prd_phi4  for
$\lambda\to 0$, i.e. infrared freedom in the present case: $\lambda(\mu\ll
m)\simeq[b_{0}\ln(m/\mu)]^{-1}$.

###  III.1 LO RGOPT

For simple illustration let us briefly recall the one-loop resultsprl_phi4 ;
prd_phi4 , thus considering the LO term in Eq. (2), including solely the first
order subtraction term, ${\cal E}_{0}=-(m^{4}/\lambda)s_{0}$, with $s_{0}$
given in Eq.(16). Performing (18), expanding to leading order $\delta^{0}$
consistently, and taking afterwards $\delta\to 1$ yields

$(4\pi)^{2}{\cal
F}^{\delta^{0}}_{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}(\frac{m}{T}).$
(22)

At this leading order Eq.(20) is satisfied exactly, so that only Eq.(19) 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}$, explicitlyTrev ; SPT3l

$\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].$
(23)

Then the exact solution of Eq. (19), using Eq.(8), is given by the self-
consistent gap equation:

$\tilde{m}^{2}=\frac{1}{2}(\frac{b_{0}}{\gamma_{0}})\>\Sigma_{R}(\tilde{m}^{2})=(4\pi)^{2}\,b_{0}\>\Sigma_{R}(\tilde{m}^{2}).$
(24)

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}},$
(25)

the free energy Eq.(22), and therefore $\tilde{m}$ in Eq.(24), 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 observedprd_phi4  between those (one-loop) RGOPT
results and the (two-loop) 2PI resummation approach in 2PI ; 2PI3loop ,
although our construction is basically very different.
To get more insight on some properties of the solution of Eq.(24), one may
conveniently use the high-temperature expansion of the relevant $J_{n}(x)$,
from Eq.(9) with $x\equiv m/T$ 444At one-loop order this approximation is
valid at the $0.1\%$ level even for $x\raisebox{-3.69899pt}{~{}\shortstack{$<$
\\\ [-0.07cm] $\sim$}}~{}1$, sufficient for our purpose since the RGOPT one-
loop solution $\tilde{m}/T$ happens to always lies in this range.. Accordingly
Eq. (24) simplifies to a quadratic equation for $x$, with a unique physical
($x>0$) solution:

$\displaystyle\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}}\simeq\pi\,\sqrt{\frac{2}{3}b_{0}\lambda}-\pi\,b_{0}\lambda+{\cal
O}(\lambda^{3/2})$ (26)

with $L_{T}\equiv\ln[\mu\,e^{\gamma_{E}}/(4\pi T)]$. Perturbatively at high
temperature, the variationally determined mass Eq.(26) has the form of a
screening mass, $\tilde{m}^{2}\sim{\cal O}(\lambda T^{2})$, but note that this
variational mass parameter is unrelated to the physical Debye mass mdeb
definition. The corresponding one-loop RGOPT pressure, from Eq.(22), 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}+{\cal
O}(\tilde{x}^{6}).$ (27)

Eqs. (24)-(27) clearly involve an all order dependence in $\lambda$, and RG-
invariance is once more manifest since from Eq.(25),
$1/(b_{0}\lambda(\mu))+L_{T}$ is $\mu$-independent. Accordingly Eqs. (26) and
(27) only depend on the single parameter $b_{0}\lambda(\mu_{0})$, where
$\mu_{0}$ is some reference scale, typically $\mu_{0}=2\pi T_{0}$.

Expanding perturbatively Eq. (27) 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}+{\cal
O}(\alpha^{5/2})$ (28)

where $\alpha\equiv b_{0}\lambda$. Note in particular that Eq.(28) contains
the nonanalytic term $\sim\lambda^{3/2}$, originating from the bosonic zero
mode resummation, but here readily obtained from RG properties. Expanding at
higher orders Eq.(26) it is easily seen that it entails nonanalytic terms
$\lambda^{(2p+1)/2},p\geq 1$ at all orders. Incidentally it is worth
mentioning that the higher orders beyond Eq.(28) obtained from Eq.(26)
correctly reproduce all orders of the $O(N$) scalar model large $N$-results
(e.g. Eq. (5.8) of phi4N ) (once including higher order ${\cal O}(x^{6})$
terms, not given in Eq.(9)), as can be checked upon identifying the correct
large-$N$ $b_{0}=1/(16\pi^{2})$ valuephi4N . 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.

Figure 2: The graphs being resummed at first nontrival RGOPT order.

###  III.2 NLO RGOPT

The NLO (two-loop) ${\cal O}(\delta^{1})$ contribution to the free energy, for
$\delta=1$, takes a rather compact form in terms of $\Sigma_{R}$ in Eq. (23):

${\cal F}^{\delta^{1}}_{0}=\frac{{\cal
E}^{\delta^{1}}_{0}}{(4\pi)^{2}}+\frac{T}{2}\hbox{$\sum$}\\!\\!\\!\\!\\!\\!\\!\int_{\bf
p}\ln(\omega_{n}^{2}+\omega_{\bf
p}^{2})|_{\mbox{finite}}-\left(\frac{2\gamma_{0}}{b_{0}}\right)\frac{m^{2}}{\lambda}\,\Sigma_{R}+\frac{\Sigma_{R}^{2}}{2\lambda},$
(29)

where the subscript on the integral term means to taking its corresponding
finite (renormalized) expression. Eq.(15) gives explicitly

${\cal
E}^{\delta^{1}}_{0}=-m^{4}\left(\frac{1}{3b_{0}\,\lambda}+\frac{s_{1}}{3}\right).$
(30)

The exact two-loop OPT and RG Eqs. (19) and (20) can be written compactly
asprd_phi4  555Note a typo in Eq.(7.2) of Ref.prd_phi4 , corrected in Eq.(31),
that does not affect any numerical results.

$f^{NLO}_{\rm
OPT}=\frac{2}{3}h\left(-s_{1}-\frac{1}{b_{0}\lambda}\right)+\frac{2}{3}S+\frac{\Sigma^{\prime}_{R}}{2m}\left(S-\frac{1}{3\lambda}\right)=0,$
(31) $f^{NLO}_{\rm
RG}=h\left[\frac{1}{6}+\left(\frac{b_{1}}{3b_{0}}-S\right)\lambda\right]+\frac{1}{2}\beta^{(2)}(\lambda)S^{2}=0,$
(32)

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)$ (33)

with, from Eq. (23),

$\Sigma^{\prime}_{R}\equiv\partial_{m^{2}}(\Sigma_{R})=\lambda(S+m^{2}S^{\prime})=\gamma_{0}\lambda\,[\ln(m^{2}/\mu^{2})-J_{2}(m/T)].$
(34)

At this stage at NLO in principle we could use three different possible
prescriptions to obtain a thermally dressed mass $\bar{m}(\lambda,T)$ as
function of the coupling, as investigated in details in prd_phi4 : Either the
OPT Eq.(31), or alternatively the (massless) RG Eq.(32), or else the full RG,
Eq.(13). The latter is not an independent equation, being a linear combination
of Eq.(31) and (32):

$f^{NLO}_{\rm full\,RG}\equiv f_{\rm RG}+2\gamma_{m}(\lambda)f_{\rm OPT}=0\;.$
(35)

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 happensrgopt1  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. (26), (27) for the first two order terms, but contain
modifications at higher ordersprd_phi4 . However, for rather large $g(2\pi
T)\lesssim 1$ reference coupling and $\mu=4\pi T$, the OPT Eq.(31) no longer
gives a real solution: in fact, the possible occurence of nonreal solutions at
higher orders is a recurrent burden of such variational approach, quite
generically expected from the nonlinear $m$ dependence if requiring to solve
$\tilde{m}(\lambda,T)$ exactly. Alternatively using Eq.(32) to determine
$\tilde{m}(\lambda,T)$ gives an unphysical solution at NLOprd_phi4 , being
driven towards the NLO UV fixed point at $\lambda=-b_{0}/b_{1}$ (an artefact
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. (13), explicitly Eq.(35) at NLO, that happens to give real
solutions at least in a much larger range of coupling values 666In RGOPT
applications to thermal QCDprdCOLD ; rgopt_NLOQCD2 , 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..

Figure 3: LO (dotted lines) and NLO (dashed, green) RGOPT pressure
$P/P_{0}(g\equiv\sqrt{\lambda/24})$ versus ${\cal O}(g^{2})$ and ${\cal
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\leq\mu\leq 4\pi T$.

The NLO pressure $P/P_{0}$ with exact $T$-dependence obtained from Eq. (35),
as function of the reference coupling $g\equiv\sqrt{\lambda(2\pi T)/24}$, is
illustrated in Fig. 3, with scale dependence $\pi T\leq\mu\leq 4\pi T$ using
the exact two-loop running Eq.(50). It is compared with LO RGOPT, Eq.(27), and
importantly also with the PT pressure, Eq.(12), respectively truncated at the
lowest ${\cal O}(\lambda)\sim g^{2}$ and next ${\cal 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 ${\cal 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\raisebox{-3.69899pt}{~{}\shortstack{$>$ \\\ [-0.07cm] $\sim$}}~{}0.6$. A
more accurate analysisprd_phi4  shows that at NLO the remnant scale dependence
reappears first at perturbative order $\lambda^{3}$. As above mentioned 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.(25) perfectly matching Eqs.(26),(27), since
the latter only depends on $b_{0}$. In contrast at NLO and beyond,
$P(\overline{m}(\lambda),T,\mu)$ inevitably has a remnant scale dependence,
basically because the subtractions in Eq.(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.(18)
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.(27), 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, conventionally often
usedTrev  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. (19) for $\tilde{\lambda}(m)$ with a resulting $P(m/T)$, and
inserting the physical Debye screening mass $m(\lambda)$ mdeb , correctly
reproducesprd_phi4  the first two terms of the standard PT pressure Eq.(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
Eq.(31) or (32), 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.SPT3l .
Note that all the relevant SPT expressions can be obtained consistently by 1)
discarding the vacuum energy subtraction ${\cal E}_{0}$ in Eq.(6); 2) taking
$a=1/2$ in Eq.(18), expanding the result to order $\delta$, and setting
$\delta\to 1$; and finally 3) calculating the variational mass gap, Eq. (19),
that gives explicitly 777Eq.(36) is what is called the tadpole prescription in
Ref.SPT3l .

$m^{2}_{\rm SPT}=\Sigma_{R},$ (36)

see Eq.(23), to be solved self-consistently for $\tilde{m}_{\rm SPT}$.
Alternatively a simpler prescription was also used in Ref.SPT3l , taking
instead the (NLO) perturbative Debye screening massmdeb :

$m^{2}_{D}=\frac{\lambda}{24}\,T^{2}\left(1-\frac{3}{\pi}(\frac{\lambda}{24})^{1/2}\right),$
(37)

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\leq\mu\leq 4\pi T$ is $\simeq
8\%$, $1.5\%$, $0.8\%$, and $0.4\%$ respectively for the ${\cal 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), if neglecting the subtraction
terms in Eq.(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 two is
not spectacular. But conversely this can largely explain why a very sizeable
scale dependence resurfaces for the NNLO SPT pressureSPT3l , 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+{\cal
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 ${\cal 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.

Figure 4: Relative scale dependence at successive orders of the running
coupling in the $\lambda\phi^{4}$ model, with $g\equiv\sqrt{\lambda/24}$ and
$\pi T\leq\mu\leq 4\pi T$ (from lower to higher values respectively).

Moreover, concerning the $\lambda\phi^{4}$ model, the peculiar sign
alternating beta-function coefficients $b_{i}$ from one- to three-loop orders,
see Eq.(45), implies that considering the running coupling alone, at three-
loop order it has a comparatively worse scale dependence than at two-loop
order, as illustrated in Fig.4. This feature tends to partly counteract the
benefits of our RG-improved construction, when comparing NLO with NNLO888We
expect the corresponding behavior in QCD to be better, since the first three
QCD beta-function coefficients have the same sign..

Applying the variational modification from Eqs.(18),(21) to the complete
three-loop free energy, sum of Eq.(6) and Eq.(10), expanded consistently now
to order $\delta^{2}$, and taking $\delta\to 1$, gives after algebra,

$\displaystyle(4\pi)^{2}{\cal F}_{0}^{\delta^{2}}$ $\displaystyle=$
$\displaystyle{\cal
E}_{0}^{\delta^{2}}-\frac{m^{4}}{8}\left(3+2L+8(\frac{\gamma_{0}}{b_{0}})^{2}(L+J_{2})\right)-\frac{T^{4}}{2}J_{0}$
(38)
$\displaystyle+\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}\left(L+J_{2}\right)+\frac{\Sigma_{R}^{2}}{128\pi^{2}\gamma^{2}_{0}\lambda}+F_{3l},$

omitting the argument $m/T$ in thermal functions $J_{i}$, with
$L\equiv\ln(\mu^{2}/m^{2})$ and $\Sigma_{R}$ defined in Eq.(23). The original
${\cal O}(\lambda^{2})$ perturbative contribution, designated by $F_{3l}$ in
Eq.(38), is already defined in Eq.(10), those terms being unaffected by
Eq.(18) at $\delta^{2}$ re-expansion order. The subtraction contributions in
Eq.(38) after the modifications from Eq.(18) read explicitly, using
$\gamma_{0}/b_{0}=1/6$:

${\cal
E}^{\delta^{2}}_{0}=-m^{4}\left(\frac{14}{81}\frac{s_{0}}{\lambda}+\frac{2}{9}s_{1}+\frac{1}{3}s_{2}\lambda+s_{3}\lambda^{2}\right).$
(39)

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

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 $\bar{m}(\lambda,T)$: namely taking the
solution of the OPT Eq.(19), or the massless RG Eq.(20), or the full RG
Eq.(13). One may further exploit the freedom to incorporate the highest order
subtraction term $s_{3}$ of Eq.(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 similarly to what happens at NLO. The behavior
of those solutions for the different prescriptions is detailed in Appendix A
for completeness.

The outcome is, in order to maximize the range of coupling and scale $\pi
T\leq\mu\leq 4\pi T$ values where real solutions are obtained, it is
appropriate to minimally neglect $s_{3}$ if using the OPT Eq.(19), and
incorporating $s_{3}\neq 0$ when using the RG Eq.(20). We mention that at NNLO
the full RG Eq.(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.(13) maximally, keeping only ${\cal 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}\neq 0$ or
taking the more involved exact Eq.(13), renders the expression even more
nonlinear in $m$, favoring the occurence of nonreal solutions. But for the
massless RG Eq.(20), provided that it has real solutions, incorporating
$s_{3}$ could be expected to give better results, since the subtraction
coefficients in Eq.(17) entering the pressure are originating directly from
the RG coefficients.

Fig. 5 illustrates as a function of the (rescaled) coupling $g$ our two
different prescriptions thus retained at NNLO: respectively applying to
Eq.(38) the OPT Eq.(19), or the RG Eq.(20) (see also Eq.(40)), with resulting
(unique) perturbatively matching exact solutions $\tilde{m}^{OPT}(\lambda,T)$
and $\tilde{m}^{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.

Figure 5: NNLO (three-loops) $\overline{m}$ from the OPT Eq.(with $s_{3}=0$)
and RG Eq (with $s_{3}\neq 0$) as function of $g\equiv\sqrt{\lambda/24}$, with
scale dependence $\pi T\leq\mu\leq 4\pi T$.

The previous exact $\tilde{m}$ solutions 999At NNLO we use the exact
expressions of all thermal integrals, in particular Eqs.(54),(57), as their
high-$T$ expressions in Eqs.(60), (61) are not very good approximations. are
then inserted within Eq.(38) to give the physical pressure. Fig. 6 illustrates
the corresponding pressures obtained from the two alternative OPT and RG mass
prescriptions. As one can see, despite the quite different $\overline{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 higher order of those a priori different prescriptions.

Figure 6: Comparison of NNLO pressure $P/P_{ideal}$ for the OPT and RG
$\overline{m}$ prescriptions (with $s_{3}=0$ and $s_{3}\neq 0$ respectively),
as function of $g\equiv\sqrt{\lambda/24}$ with $\pi T\leq\mu\leq 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}$.

Finally we illustrate our main results at successive LO, NLO, and NNLO in Fig.
7, compared with the PT pressure, Eq.(12), and SPT pressure at successive
orders. The NNLO PT, with successive terms up to ${\cal O}(\lambda^{2})$, has
a substantially larger remnant scale dependence than the ${\cal
O}(\lambda^{3/2})$ PT in Fig. 3. Concerning the SPT, at NNLO we only show the
results from using the mass gap Eq.(36), following the very same prescription
as in SPT3l , namely using also Eq.(36) at NNLO. (N.B. using the perturbative
screeening mass instead, Eq.(37), 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 moderate
remnant scale dependence at NNLO, almost invisible at the figure scale until
relatively large $g\simeq 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 counter-effect 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.

Figure 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\leq\mu\leq 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}$.

##  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.(20), thus
more directly embedding RG properties, appears to give a slightly better
remnant scale dependence than the more traditional variational mass
prescription from Eq.(19) optimizing the pressure.

The RGOPT has been applied recently at NLO in thermal QCDrgopt_NLOQCD1 ;
rgopt_NLOQCD2 , 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 pseudo-transition, 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.

##  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.(20), straightforward to obtain after algebra from Eq.(38), that
takes the rather compact form:

$\displaystyle f^{NNLO}_{\rm RG}$ $\displaystyle=$
$\displaystyle-\frac{m^{4}}{2}\left(1-6\frac{\gamma_{0}}{b_{0}}+8(\frac{\gamma_{0}}{b_{0}})^{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}$
(40)
$\displaystyle+\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}{81}\frac{s_{0}}{\lambda^{2}}+\frac{s_{2}}{3}+2s_{3}\lambda\right)\right\\}=0,$

where

$\displaystyle\displaystyle
G_{3l}=-\frac{1}{24}(\frac{\lambda}{16\pi^{2}})^{2}$
$\displaystyle\Big{[}m^{4}\left(15L^{2}+34L+\frac{41}{2}+6(L+1)J_{2}\right)$
(41)
$\displaystyle-m^{2}T^{2}J_{1}\left(24L+28+6J_{2}\right)+9T^{4}J^{2}_{1}\Big{]},$
$\beta^{(3)}(\lambda)=\lambda^{2}(b_{0}+b_{1}\lambda+b_{2}\lambda^{2}),$ (42)

with $F_{3l}$ given in Eq.(10) and $\Sigma_{R}$ in Eq.(23). The OPT Eq.(19) at
NNLO is similarly easily obtained from Eq.(38), but it gives a somewhat more
lengthy expression that we thus refrain to give explictly.

Figure 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$.

Fig. 8 illustrates, for $s_{3}=0$ (top) and $s_{3}\neq 0$ (bottom), the
behavior of these exact RG Eq.(40) and OPT (19) 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$), 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 (bottom). 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 behaviour is obtained but the disappearance of real
perturbative-matching RG solution is delayed to substantially higher $g\gtrsim
1$ values for any $\mu\leq 4\pi T$: accordingly for $0<g\leq 1$ and $\pi
T\leq\mu\leq 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 $\overline{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.(13) has no real solutions for $0<g\leq 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}+\cdots$
(43)

and

$\gamma_{m}(\lambda)\equiv\frac{d\ln
m}{d\ln\mu}=\gamma_{0}\lambda+\gamma_{1}\lambda^{2}+\gamma_{2}\lambda^{3}+\cdots$
(44)

where up to three-loop order RGphi4loop

$(4\pi)^{2}b_{0}=3;\;\;(4\pi)^{4}\,b_{1}=-17/3;\;(4\pi)^{6}\,b_{2}=\frac{3915+2592\zeta(3)}{216},$
(45)
$(4\pi)^{2}\gamma_{0}=1/2;\;\;(4\pi)^{4}\,\gamma_{1}=-5/12;\;\;(4\pi)^{6}\,\gamma_{2}=\frac{7}{4}.$
(46)

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 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}+{\cal
O}(\lambda^{3})$ (47)
$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}+{\cal
O}(\lambda^{3})$ (48)

The necessary additional vacuum energy counterterms up to three loop order
readvacanom_kastening  in our normalization conventions:

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

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

$\lambda^{\rm(2-loop)}(\mu)=\frac{\lambda(\mu_{0})}{f_{W}\left(\lambda(\mu_{0}),\ln\frac{\mu}{\mu_{0}}\right)}\;\;,$
(50)

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^{-1}_{W}}{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^{-[1+\frac{b_{0}}{b_{1}\lambda}(1-b_{0}\lambda
L)]}\right]\right\\}\;\;,$ (51)

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.(50), (51) give not much
visible differences with a simpler perturbatively truncated expansion at order
$\lambda^{3}$:

$\lambda^{-1}(\mu)\simeq\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}+{\cal
O}(\lambda^{4})\;.$ (52)

At 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 not too large
coupling values there is not much visible differences with a more common
perturbatively truncated running coupling:

$\displaystyle\lambda^{\rm(3-loop)}(\mu)=\frac{\lambda(\mu_{0})}{f^{3L}},$
$\displaystyle f^{3L}$ $\displaystyle=$ $\displaystyle 1-b_{0}L_{\mu}\lambda-
b_{1}L_{\mu}\lambda^{2}-(b_{0}\frac{b_{1}}{2}L_{\mu}^{2}+b_{2}L_{\mu})\lambda^{3}$
$\displaystyle+$
$\displaystyle\frac{L_{\mu}}{6b_{0}^{2}}\left(6b_{1}(b_{1}^{2}-2b_{0}b_{2})-3b_{0}^{2}(b_{1}^{2}+2b_{0}b_{2})L_{\mu}-2b_{0}^{4}b_{1}L_{\mu}^{2}\right)\lambda^{4}.$
(53)

##  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
scal3l , and entering Eq.(10):

$K_{2}(\frac{m}{T})=-\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),$
(54)

where

$\displaystyle f_{2}(E,k)$
$\displaystyle=\left(\frac{E^{2}-M_{k}^{2}}{E^{2}-k^{2}}\right)^{\frac{1}{2}}\ln\frac{\left(E^{2}-k^{2}\right)^{\frac{1}{2}}+\left(E^{2}-M_{k}^{2}\right)^{\frac{1}{2}}}{\left(E^{2}-k^{2}\right)^{\frac{1}{2}}-\left(E^{2}-M_{k}^{2}\right)^{\frac{1}{2}}}\
\ \ \ ,\ \ k^{2}<E^{2}-4m^{2}$ (55)
$\displaystyle=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}$
$\displaystyle=\left(\frac{M_{k}^{2}-E^{2}}{k^{2}-E^{2}}\right)^{\frac{1}{2}}\ln\frac{\left(M_{k}^{2}-E^{2}\right)^{\frac{1}{2}}+\left(k^{2}-E^{2}\right)^{\frac{1}{2}}}{\left(M_{k}^{2}-E^{2}\right)^{\frac{1}{2}}-\left(k^{2}-E^{2}\right)^{\frac{1}{2}}}\
\ \ \ ,\ \ E^{2}<k^{2}$

and

$\displaystyle M_{k}^{2}$ $\displaystyle=$ $\displaystyle
4m^{2}+k^{2},\;\;k\equiv|{\bf p}+{\bf q}|,$ $\displaystyle E_{\sigma}(p,q)$
$\displaystyle=$ $\displaystyle\sqrt{p^{2}+m^{2}}+\sigma\sqrt{q^{2}+m^{2}}.$
(56) $\displaystyle K_{3}(\frac{m}{T})=\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}}$
$\displaystyle\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)\\},$
(57)

where

$\displaystyle f_{3}(E,p)$
$\displaystyle=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}$ (58)
$\displaystyle=p\ln\frac{\big{|}E^{2}-m^{2}-p^{2}\big{|}}{m^{2}}+\left(E^{2}-m^{2}\right)^{\frac{1}{2}}\ln\frac{(E^{2}-m^{2})^{\frac{1}{2}}+p}{\big{|}(E^{2}-m^{2})^{\frac{1}{2}}-p\big{|}}\
\ \ ,\ E^{2}>m^{2}$

and

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

In the limit $x\equiv m/T\rightarrow 0$ Eqs.(54), (57) can be expressed
analytically as:

$K_{2}(x)\simeq\frac{(4\pi)^{4}}{72}\left(\ln
x+\frac{1}{2}+\frac{\zeta^{\prime}(-1)}{\zeta(-1)}\right)-372.65\,x\left(\ln
x+1.4658\right),$ (60)
$K_{3}(x)\simeq\frac{(4\pi)^{4}}{48}\left(-\frac{7}{15}+\frac{\zeta^{\prime}(-1)}{\zeta(-1)}-\frac{\zeta^{\prime}(-3)}{\zeta(-3)}\right)+1600.0\,x\left(\ln
x+1.3045\right).$ (61)

For the numerical results we rather use the exact expressions Eqs.(54), (57).

## References

  * (1) J. P. Blaizot, E. Iancu and A. Rebhan, In *Hwa, R.C. (ed.) et al.: Quark gluon plasma* 60-122 [hep-ph/0303185]; U. Kraemmer and A. Rebhan, Rept. Prog. Phys. 67, 351 (2004) [hep-ph/0310337].
  * (2) J. I. Kapusta and C. Gale, “Finite-temperature field theory: Principles and applications” (Cambridge University Press, 2006).
  * (3) M. Laine and A. Vuorinen, Lect. Notes Phys. 925, 1 (2016). [arXiv:1701.01554 [hep-ph]].
  * (4) Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz and K. K. Szabo, Nature 443, 675-678 (2006) [arXiv:hep-lat/0611014 [hep-lat]]; Y. Aoki, Z. Fodor, S. D. Katz and K. K. Szabo, Phys. Lett. B 643, 46-54 (2006) [arXiv:hep-lat/0609068 [hep-lat]]; M. Cheng et al, et al. Phys. Rev. D 74, 054507 (2006) [arXiv:hep-lat/0608013 [hep-lat]].
  * (5) P. de Forcrand, PoS LAT 2009, 010 (2009); [arXiv:1005.0539 [hep-lat]]; G. Aarts, J. Phys. Conf. Ser. 706, 022004 (2016). [arXiv:1512.05145 [hep-lat]].
  * (6) R. R. Parwani, Phys. Rev. D 45, 4695 (1992) [erratum: Phys. Rev. D 48, 5965 (1993)] [arXiv:hep-ph/9204216 [hep-ph]].
  * (7) F. Karsch, A. Patkos and P. Petreczky, Phys. Lett. B 401, 69-73 (1997) [arXiv:hep-ph/9702376 [hep-ph]].
  * (8) J. P. Blaizot, A. Ipp and N. Wschebor, Nucl. Phys. A 849, 165-181 (2011) [arXiv:1007.0991 [hep-ph]].
  * (9) J. P. Blaizot, E. Iancu and A. Rebhan, Phys. Rev. D 63, 065003 (2001) [arXiv:hep-ph/0005003 [hep-ph]].
  * (10) J. Berges, S. Borsanyi, U. Reinosa and J. Serreau, Phys. Rev. D 71, 105004 (2005) [arXiv:hep-ph/0409123 [hep-ph]]; J. Berges, S. Borsanyi, U. Reinosa and J. Serreau, Annals Phys. 320, 344-398 (2005) [arXiv:hep-ph/0503240 [hep-ph]].
  * (11) M. E. Carrington, B. A. Meggison and D. Pickering, Phys. Rev. D 94, no.2, 025018 (2016) [arXiv:1603.02085 [hep-th]].
  * (12) J. P. Blaizot and N. Wschebor, Phys. Lett. B 741, 310-315 (2015) [arXiv:1409.4795 [hep-ph]].
  * (13) J. O. Andersen, E. Braaten and M. Strickland, Phys. Rev. D 63, 105008 (2001) [arXiv:hep-ph/0007159 [hep-ph]].
  * (14) J. O. Andersen and M. Strickland, Phys. Rev. D 64, 105012 (2001) [arXiv:hep-ph/0105214 [hep-ph]]; J. O. Andersen and M. Strickland, Annals Phys. 317, 281-353 (2005) [arXiv:hep-ph/0404164 [hep-ph]].
  * (15) R. L. S. Farias, G. Krein and R. O. Ramos, Phys. Rev. D 78, 065046 (2008) [arXiv:0809.1449 [hep-ph]].
  * (16) J. O. Andersen and L. Kyllingstad, Phys. Rev. D 78, 076008 (2008) [arXiv:0805.4478 [hep-ph]].
  * (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. B85, 421 (1979); R.P. Feynman and H. Kleinert, Phys. Rev. A34, 5080 (1986); A. Duncan and M. Moshe, Phys. Lett. B 215, 352-358 (1988) H.F. Jones and M. Moshe, Phys. Lett. B234, 492 (1990); A. Neveu, Nucl. Phys. B, Proc. Suppl. B18, 242 (1991); V. Yukalov, J. Math. Phys (N.Y) 32, 1235 (1991); S. Gandhi, H.F. Jones and M. Pinto, Nucl. Phys. B359, 429 (1991); C. M. Bender et al., Phys. Rev. D45, 1248 (1992); H. Yamada, Z. Phys. C59, 67 (1993); A.N. Sisakian, I.L. Solovtsov and O.P. Solovtsova, Phys. Lett. B321, 381 (1994); H. Kleinert, Phys. Rev. D57, 2264 (1998); Phys. Lett. B434, 74 (1998).
  * (18) R. Seznec and J. Zinn-Justin, J. Math. Phys. 20, 1398 (1979); J.C. Le Guillou and J. Zinn-Justin, Ann. Phys. 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, 1827 (1992).
  * (21) J. O. Andersen, E. Braaten and M. Strickland, Phys. Rev. Lett. 83, 2139 (1999) [hep-ph/9902327]; J. O. Andersen, E. Braaten and M. Strickland, Phys. Rev. D 61, 074016 (2000) [arXiv:hep-ph/9908323 [hep-ph]].
  * (22) J. O. Andersen, M. Strickland and N. Su, Phys. Rev. Lett. 104, 122003 (2010) [arXiv:0911.0676 [hep-ph]]; JHEP 08, 113 (2010) [arXiv:1005.1603 [hep-ph]]; J. O. Andersen, L. E. Leganger, M. Strickland and N. Su, JHEP 08, 053 (2011) [arXiv:1103.2528 [hep-ph]].
  * (23) S. Mogliacci, J. O. Andersen, M. Strickland, N. Su and A. Vuorinen, JHEP 12, 055 (2013) [arXiv:1307.8098 [hep-ph]]; N. Haque, J. O. Andersen, M. G. Mustafa, M. Strickland and N. Su, Phys. Rev. D 89, no.6, 061701 (2014) [arXiv:1309.3968 [hep-ph]]; N. Haque, A. Bandyopadhyay, J. O. Andersen, M. G. Mustafa, M. Strickland and N. Su, JHEP 05, 027 (2014) [arXiv:1402.6907 [hep-ph]].
  * (24) J. L. Kneur and A. Neveu, Phys. Rev. D 81, 125012 (2010) [arXiv:1004.4834 [hep-th]].
  * (25) J.-L. Kneur and A. Neveu, Phys. Rev. D 88, no.7, 074025 (2013) [arXiv:1305.6910 [hep-ph]].
  * (26) J.-L. Kneur and A. Neveu, Phys. Rev. D 92, no.7, 074027 (2015) [arXiv:1506.07506 [hep-ph]].
  * (27) J. L. Kneur and A. Neveu, Phys. Rev. D 101, no.7, 074009 (2020) [arXiv:2001.11670 [hep-ph]].
  * (28) J. L. Kneur and M. B. Pinto, Phys. Rev. D 92, 116008 (2015) [arXiv:1508.02610 [hep-ph]].
  * (29) J. L. Kneur and M. B. Pinto, Phys. Rev. Lett. 116, 031601 (2016) [arXiv:1507.03508 [hep-ph]].
  * (30) G. N. Ferrari, J. L. Kneur, M. B. Pinto and R. O. Ramos, Phys. Rev. D 96, 116009 (2017) [arXiv:1709.03457 [hep-ph]].
  * (31) J. L. Kneur, M. B. Pinto and T. E. Restrepo, Phys. Rev. D 100, 114006 (2019) [arXiv:1908.08363 [hep-ph]].
  * (32) J. L. Kneur, M. B. Pinto and T. E. Restrepo, [arXiv:2101.02124 [hep-ph]] (to be published in Phys.Rev. D letter.)
  * (33) J. L. Kneur, M. B. Pinto and T. E. Restrepo, [arXiv:2101.08240 [hep-ph]] (to be published in Phys.Rev. D.)
  * (34) J. O. Andersen, E. Braaten and M. Strickland, Phys. Rev. D 62, 045004 (2000) [arXiv:hep-ph/0002048 [hep-ph]].
  * (35) J. Frenkel, A. V. Saa and J. C. Taylor, Phys. Rev. D 46, 3670-3673 (1992); P. B. Arnold and C. X. Zhai, Phys. Rev. D 50, 7603-7623 (1994) [arXiv:hep-ph/9408276 [hep-ph]]; R. Parwani and H. Singh, Phys. Rev. D 51, 4518-4524 (1995) [arXiv:hep-th/9411065 [hep-th]].
  * (36) E. Braaten and A. Nieto, Phys. Rev. D 51, 6990-7006 (1995) [arXiv:hep-ph/9501375 [hep-ph]].
  * (37) B. M. Kastening, Phys. Rev. D 54, 3965-3975 (1996) [arXiv:hep-ph/9604311 [hep-ph]]; Phys. Rev. D 57, 3567-3578 (1998) [arXiv:hep-ph/9710346 [hep-ph]].
  * (38) I. T. Drummond, R. R. Horgan, P. V. Landshoff and A. Rebhan, Nucl. Phys. B 524, 579-600 (1998) [arXiv:hep-ph/9708426 [hep-ph]].
  * (39) E. Brezin, J. C. Le Guillou and J. Zinn-Justin, Phys. Rev. D 9, 1121 (1974).