# Chiral Condensate and Spectral Density at full five-loop and partial six-loop orders of Renormalization Group Optimized Perturbation

Jean-Loïc Kneur, Andre Neveu

► To cite this version:

Jean-Loïc Kneur, Andre Neveu. Chiral Condensate and Spectral Density at full five-loop and partial six-loop orders of Renormalization Group Optimized Perturbation. Physical Review D, 2020, pp.074009. 10.1103/PhysRevD.101.074009. hal-02464734

HAL Id: hal-02464734

https://hal.science/hal-02464734

Submitted on 4 Jun 2021

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
# Chiral condensate and spectral density at full five-loop and partial six-loop orders of renormalization group optimized perturbation theory

Jean-Loïc Kneur and André Neveu

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

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

(Received 6 February 2020; accepted 9 March 2020; published 9 April 2020)

We reconsider our former determination of the chiral quark condensate $\langle \bar{q}q \rangle$ from the related QCD spectral density of the Euclidean Dirac operator, using our renormalization group optimized perturbation (RGOPT) approach. Thanks to the recently available complete five-loop QCD RG coefficients, and some other related four-loop results, we can extend our calculations exactly to $N^4LO$ (five loops) RGOPT, and partially to $N^5LO$ (six loops), the latter within a well-defined approximation accounting for all six-loop contents exactly predictable from five-loop RG properties. The RGOPT results overall show a very good stability and convergence, giving primarily the RG-invariant (RGI) condensate, $\langle \bar{q}q \rangle_{\text{RGI}}^{1/3}(n_f=0) = -(0.840^{+0.020}_{-0.016})\bar{\Lambda}_0$, $\langle \bar{q}q \rangle_{\text{RGI}}^{1/3}(n_f=2) = -(0.781^{+0.019}_{-0.009})\bar{\Lambda}_2$, $\langle \bar{q}q \rangle_{\text{RGI}}^{1/3}(n_f=3) = -(0.751^{+0.019}_{-0.010})\bar{\Lambda}_3$, where $\bar{\Lambda}_{n_f}$ is the basic QCD scale in the $\overline{\text{MS}}$ scheme for $n_f$ quark flavors, and the range spanned is our rather conservative estimated theoretical error. This leads e.g., to $\langle \bar{q}q \rangle_{n_f=3}^{1/3}(2 \text{ GeV}) = -(273^{+7}_{-4} \pm 13) \text{ MeV}$, using the latest $\bar{\Lambda}_3$ values giving the second uncertainties. We compare our results with some other recent determinations. As a by-product of our analysis we also provide complete five-loop and partial six-loop expressions of the perturbative QCD spectral density, that may be useful for other purposes.

DOI: 10.1103/PhysRevD.101.074009

## I. INTRODUCTION

The chiral quark condensate $\langle \bar{q}q \rangle$ is a main order parameter of spontaneous chiral symmetry breaking, $SU(n_f)_L \times SU(n_f)_R \to SU(n_f)_V$ for $n_f$ massless quarks. It is an intrinsically nonperturbative quantity, indeed vanishing at any finite order of ordinary perturbative QCD in the chiral limit. For nonvanishing quark masses, the famous Gell-Mann-Oakes-Renner (GMOR) relation [1], e.g., for the two lightest flavors,

$$F_{\pi}^{2} m_{\pi}^{2} = -(m_{u} + m_{d}) \langle \bar{u}u \rangle + \mathcal{O}(m_{\bar{q}}^{2}) \quad (1.1)$$

relates the condensate with the pion mass $m_\pi$ and decay constant $F_\pi$ together with the (current) quark masses. At present the light quark masses $m_{u,d,s}$ determined from lattice simulations (see [2] for a recent review) give an indirect determination of the condensate from using (1.1). Phenomenological values of the condensate can also be extracted [3,4] indirectly from data using spectral QCD sum rule methods [5]. However, the GMOR relation (1.1)

entails explicit chiral symmetry breaking from quark masses, and is valid up to higher-order terms $\mathcal{O}(m_\bar{q}^2)$. Thus more direct “first principles” determinations are always desirable to disentangle quark current mass effects for a better understanding of the dynamical chiral symmetry breaking mechanism at work in QCD. Analytical determinations have been derived in various models and approximations, starting early with the Nambu and Jona-Lasinio model [6,7]. There is also a long history of determinations based on Schwinger-Dyson equations and related approaches [8–11] typically. Lattice calculations have also determined the quark condensate by different approaches [12], in particular by computing the spectral density of the Dirac operator [13–15] directly related to the quark condensate via the Banks-Casher relation [16–18]. Although some of the lattice determinations are very precise, those always rely on extra assumptions and modelization to extrapolate to the chiral limit [19], using mainly chiral perturbation theory [20]. Moreover the convergence properties of chiral perturbation [21] for $n_f = 3$ are not as good as for $n_f = 2$, and different recent lattice simulations still show rather important discrepancies [2]. Also, within an extended chiral perturbation framework, it has been found significant suppression of the three-flavor case with respect to the two-flavor case [22], which may be attributed to the relatively large explicit chiral symmetry breaking from the strange quark mass.

<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>
Our renormalization group optimized perturbation (RGOPT) approach [23–25] provides analytic sequences of (variational) nonperturbative approximations having a nontrivial chiral limit. As such it provides in particular an alternative independent determination of the chiral condensate [26]. More generally the RGOPT method has also been explored so far in various models, in particular to improve the resummation properties of thermal perturbative expansions for thermodynamical quantities at finite temperatures [27,28], and for QCD at finite densities [29]. In the present work we iterate on our previous three- and four-loop RGOPT determination [26] of the condensate in vacuum from the related spectral density, by going at the complete five-loop and partial six-loop level of our approximation.

In Sec. II we shortly recall the well-known connection of the condensate with the spectral density of the Dirac operator through the Banks-Casher relation. Also for completeness, in Sec. III we shortly review our RGOPT variational construction of nonperturbative approximations, and its adaptation to the evaluation of the spectral density, as already detailed in Ref. [26]. In Sec. IV we derive the standard perturbative quark condensate and related perturbative spectral density, exactly up to five-loop order and partially up to six-loop order in a well-defined approximation, thanks most notably to the recently available five-loop RG coefficients [30,31], in particular the crucially relevant vacuum anomalous dimension [32]. The perturbative spectral density for an arbitrary number of quark flavors can also be useful for other purposes irrespective of our variational approach, most typically for perturbative matching of lattice simulation results. Section V gives our detailed numerical analysis and the RGOPT condensate results order by order up to five and (approximate) six loops, discussing also different approximation variants in order to estimate the theoretical uncertainties of our predictions. In Sec. VI we compare with other recent determinations, mainly from lattice simulations. Finally Sec. VII presents a summary and conclusions, and the Appendix completes various relevant expressions.

## II. SPECTRAL DENSITY AND THE QUARK CONDENSATE

For a more detailed review of the connection of the density of eigenvalues $\rho(\lambda)$ of the Dirac operator with the chiral condensate $\langle \bar{q}q \rangle$ through the Banks-Casher relation [16], we refer to our previous four-loop analysis [26] and to former works and reviews (see e.g., [17]). The link between the spectral density and the condensate appearing in the operator product expansion has been carefully discussed in [10]. In short, in the infinite volume limit the spectrum of the Euclidean Dirac operator becomes dense, and using the formal definition of the quark condensate together with the properties of the eigenvalues of the Dirac operator leads to the relation

$$ \langle \bar{q}q \rangle (m) = -2m \int_0^\infty d\lambda \frac{\rho(\lambda)}{\lambda^2 + m^2}. \quad (2.1) $$

Equation (2.1) essentially expresses that the two-point quark correlator has a spectral representation as a function of $m$. The Banks-Casher relation is the chiral symmetric limit of Eq. (2.1), that gives the chiral condensate as

$$ \lim_{m \to 0} \langle \bar{q}q \rangle = -\pi \rho(0), \quad (2.2) $$

if the spectral density at the origin can be determined. Note also that for nonzero fermion mass $m$, the spectral density is thus determined by the discontinuity of $\langle \bar{q}q \rangle (m)$ across the imaginary axis:

$$ \rho(\lambda) = -\frac{1}{2\pi} [\langle \bar{q}q \rangle (i\lambda + \epsilon) - \langle \bar{q}q \rangle (i\lambda - \epsilon)]|_{\epsilon \to 0}. \quad (2.3) $$

For nonvanishing quark mass $m$, $\langle \bar{q}q \rangle$ has a nontrivial perturbative series expansion, $\sim m^3 f[\ln(m^2/\mu^2)]$, and its discontinuities are simply given by those coming from the perturbative logarithmic mass dependence. Therefore the above relation (2.3) also allows us to calculate the corresponding perturbative spectral density. However, the $\lambda \to 0$ limit, relevant for the true chiral condensate, trivially leads to a vanishing result, since perturbatively $\rho(\lambda) \sim \lambda^3$. But as we recall below a crucial feature of the variational RGOPT method is to circumvent this, giving a nontrivial result for $\lambda \to 0$.

## III. RG OPTIMIZED PERTURBATION

### A. Optimized perturbation (OPT) and RGOPT construction

The RGOPT is basically a variational approach made compatible with RG properties. The starting point is to deform the standard QCD Lagrangian by introducing a variational (quark) mass term partly treated as an interaction term. One can most conveniently organize this systematically at arbitrary perturbative orders, by introducing a new expansion parameter $0 < \delta < 1$ interpolating between the (massive) free Lagrangian $\mathcal{L}_{\text{free}}$ and the original (massless) Lagrangian $\mathcal{L}_{\text{int}}$ respectively. This amounts first to the prescription

$$ m_q \to m(1-\delta)^a, \quad g \to \delta g, \quad (3.1) $$

within some given (renormalized) perturbative expansion of a physical quantity $P(m, g)$ (here $g \equiv 4\pi\alpha_S$ for QCD). In Eq. (3.1) we introduce for more generality an extra exponent $a$, that plays a crucial role in our approach, as we recall below. Next the resulting expression is expanded in powers of $\delta$ at order $k$, the so-called $\delta$ expansion [33], and afterwards $\delta \to 1$ is taken to recover the original massless theory. This leaves a remnant $m$ dependence at any finite $k$ order: Since at infinite $k$ order there is in principle no dependence on $m$, a finite-order approximation
can be obtained through an optimization (OPT) prescription, i.e., a minimization of the dependence on $m$,

$$ \frac{\partial}{\partial m} P^{(k)}(m, g, \delta = 1)|_{m \equiv \tilde{m}} = 0, \qquad (3.2) $$

determining a nontrivial dressed mass $\tilde{m}(g)$. The prescription is consistent with renormalizability [34–36] and gauge invariance, and (3.2) realizes dimensional transmutation, in contrast with the original mass vanishing in the chiral limit. In simpler one-dimensional models the procedure is a particular case of “order-dependent mapping” [37], and was shown to converge exponentially fast for the oscillator energy levels [38].

Now in most previous OPT applications, the simple (linear) value $a = 1$ was used in Eq. (3.1) for the $\delta$ expansion mainly for simplicity. In contrast we combine [23–25] the OPT Eq. (3.2) with RG properties, by requiring the ($\delta$-modified) expansion to satisfy, in addition to Eq. (3.2), a perturbative RG equation:

$$ \mu \frac{d}{d\mu} (P^{(k)}(m, g, \delta = 1)) = 0, \qquad (3.3) $$

where the (homogeneous) RG operator acting on a physical quantity is defined as<fnref n="1" />

$$ \mu \frac{d}{d\mu} = \mu \frac{\partial}{\partial\mu} + \beta(g) \frac{\partial}{\partial g} - \gamma_m(g) m \frac{\partial}{\partial m}. \qquad (3.4) $$

Note that once combined with Eq. (3.2), the RG equation takes a reduced massless form:

$$ \left[ \mu \frac{\partial}{\partial \mu} + \beta(g) \frac{\partial}{\partial g} \right] P^{(k)}(m, g, \delta = 1) = 0. \qquad (3.5) $$

Then a crucial observation is that after performing (3.1), perturbative RG invariance is generally lost, so that Eq. (3.5) gives a nontrivial additional constraint,<fnref n="2" /> but RG invariance can only be restored for a unique value of the exponent $a$ fully determined by the (scheme-independent) first-order RG coefficients [24,25]:

$$ a \equiv \gamma_0 / (2b_0). \qquad (3.6) $$

Therefore Eqs. (3.5), (3.6) and (3.2) together completely fix optimized $m \equiv \tilde{m}$ and $g \equiv \tilde{g}$ values. Moreover the prescription with (3.6) drastically improves the convergence properties [25].

Another known issue of standard OPT is that Eq. (3.2) alone generally gives more and more solutions as one proceeds to higher orders, with some being complex. Thus

it may be difficult to select the right solutions, and unphysical (nonreal) ones are a burden. In contrast, the additional constraint (3.6) guarantees that at arbitrary $\delta$ orders at least one of both the RG and OPT solutions $\tilde{g}(m)$ continuously matches the standard perturbative RG behavior for $g \to 0$ [i.e., asymptotic freedom (AF) for QCD],

$$ \tilde{g}(\mu \gg \tilde{m}) \sim \left( 2b_0 \ln \frac{\mu}{\tilde{m}} \right)^{-1} + \mathcal{O}\left( \left( \ln \frac{\mu}{\tilde{m}} \right)^{-2} \right), \qquad (3.7) $$

and these AF-matching solutions are often unique at a given order for both the RG and OPT equations. However, (3.6) does not guarantee in general that the compelling AF-matching solution remains real valued for all physically relevant ranges. Actually the occurrence of complex solutions is merely a consequence of solving exactly the (polynomial) Eqs. (3.2) and (3.5), but since those equations are derived from a perturbative expansion originally, they cannot be considered truly exact. Thus in practice one can often recover real solutions by considering a more approximate (perturbatively consistent) RG equation or solution (see e.g., [25,29]).

## B. RGOPT for the spectral density

As shortly reviewed above in Sec. II, using the spectral density with the Banks-Casher relation (2.2) gives direct access to the QCD condensate in the chiral limit. Therefore the spectral density constitutes a particularly suitable ansatz to apply our variational approach (see [26] for more discussions). The RG equation relevant for $\rho(\lambda, a_s \equiv \alpha_s/\pi)$ was derived in [26] and is completely analogous to the standard RG equation, but with the mass replaced by the spectral parameter,

$$ \left[ \mu \frac{\partial}{\partial \mu} + \beta(a_s) \frac{\partial}{\partial a_s} - \gamma_m(a_s) \lambda \frac{\partial}{\partial \lambda} - \gamma_m(a_s) \right] \rho(\lambda, a_s) = 0. \qquad (3.8) $$

One can next proceed to the modification of the resulting perturbative series $\rho(\lambda, a_s)$ as implied by the $\delta$ expansion, now, from Eq. (3.8) clearly applied not on the original mass but on the spectral value<fnref n="3" /> $\lambda$:

$$ \lambda \to \lambda(1 - \delta)^a \qquad a_s \to \delta a_s. \qquad (3.9) $$

Consequently the mass optimization on $\langle \bar{q}q \rangle$ thus translates into an optimization of the spectral density with respect to $\lambda$,

$$ \frac{\partial \rho^{(k)}(\lambda, a_s)}{\partial \lambda} = 0, \qquad (3.10) $$

at successive $\delta^k$ order (see [26] for more details).

<footnote n="1">¹Our normalization is $g \equiv 4\pi\alpha_s$, $\beta(g) \equiv dg/d \ln \mu$, $\gamma_m(g) \equiv -d \ln m/d \ln \mu$; see Appendix for relations to [30,31].</footnote>

<footnote n="2">²A connection of the exponent $a$ with RG anomalous dimensions/critical exponents had also been established previously in the $D=3 \Phi^4$ model for the Bose-Einstein condensate critical temperature shift by two independent OPT approaches [39,40].</footnote>

<footnote n="3">³We simplify notations with $\lambda \equiv |\lambda|$ since it is necessarily positive.</footnote>
Finally, as one last subtlety, note that the interpolation exponent $a$ in Eq. (3.6) is universal insofar as the original expansion to be modified is itself (perturbatively) RG invariant. Now since $m\langle\bar{q}q\rangle$ is the RG-invariant quantity, rather than $\langle\bar{q}q\rangle$, when performing the perturbative modification implied by (3.9) on the spectral density, it is easily derived that the consistent value to be used is rather

$$a = \frac{4}{3} \left( \frac{\gamma_0}{2b_0} \right), \qquad (3.11)$$

which also maintains the occurrence of essentially unique AF-matching solutions with a behavior similar to (3.7) (with $m \to \lambda$ understood).

# IV. PERTURBATIVE QUARK CONDENSATE AND SPECTRAL DENSITY

## A. Perturbative quark condensate

The perturbative expansion of the QCD quark condensate for a nonzero quark mass can be calculated systematically from the directly related vacuum-energy graphs. A few representative Feynman graph contributions at successive orders are illustrated (up to three-loop order only) in Fig. 1. (There are evidently some more three-loop contributions, not shown here.) Note that the one-loop order is $O(g^0) = O(1)$. The perturbative series for the renormalized quantity $m\langle\bar{q}q\rangle$ up to six-loop order reads formally,

$$m\langle\bar{q}q\rangle(m, a_s) = \left(\frac{3}{2\pi^2}\right) m^4 \left\{ \begin{aligned} & \frac{1}{2} - L_m + 4a_s \left(L_m^2 - \frac{5}{6}L_m + \frac{5}{12}\right) \\ & + a_s^2 \sum_{i=0}^3 c_{3i} L_m^{3-i} + a_s^3 \sum_{i=0}^4 c_{4i} L_m^{4-i} \\ & + a_s^4 \sum_{i=0}^5 c_{5i} L_m^{5-i} + a_s^5 \sum_{i=0}^6 c_{6i} L_m^{6-i} \end{aligned} \right\}, \quad (4.1)$$

where $L_m \equiv \ln(m/\mu)$, $m \equiv m(\mu)$ and $a_s \equiv \alpha_s(\mu)/\pi$ in the $\overline{\text{MS}}$ scheme with renormalization scale $\mu$. The two-loop contributions were calculated in the $\overline{\text{MS}}$ scheme long ago, first in [41] (see also [36]). At higher $k$-loop orders ($k \ge 3$) we have formally defined the coefficients as $c_{ki}$ for convenience, with their explicit expressions given below and in the Appendix. Before detailing these expressions, we recall some rather well-known but important features related to RG properties. First, note that the calculation of the graphs in Fig. 1 still contains divergent terms, not canceled by mass and coupling renormalization (as is clear already from the very first one-loop graph). Those divergences need an additive renormalization; in other words $m\langle\bar{q}q\rangle$ has its own anomalous dimension directly related to the (quark part of) vacuum-energy anomalous dimension. This also implies that the finite expression (4.1) is not separately RG invariant: More precisely the perturbative RG invariance is expressed in our normalization as

<figure><img src="image_3.png" /><figcaption>FIG. 1. Samples of standard perturbative QCD contributions to the chiral condensate up to three-loop order. The cross denotes a mass insertion.</figcaption></figure>

$$\mu \frac{d}{d\mu} (m\langle\bar{q}q\rangle(m, a_s)) + 4m^4\Gamma^0(a_s) \equiv 0. \quad (4.2)$$

where the first term is the (homogeneous) RG operator given in Eq. (3.4) and $\Gamma^0$ is the vacuum-energy anomalous dimension [41,42], remarkably recently evaluated fully analytically to five loops by the authors of Ref. [32] [see more details in Eq. (A2) in the Appendix]. Therefore note that the RG consistency expressed by requiring Eq. (4.2) to hold perturbatively order by order allows us to determine all the logarithmic $(L_m)^p$ coefficients $c_{ki}$, with $k \ge i + 2$ at perturbative orders $k$ from lower ($< k$) order coefficients and RG $\beta(a_s)$ and $\gamma_m(a_s)$ functions up to order $k - 2$ and $k - 1$ respectively. In addition the knowledge of $\Gamma^0(a_s)$ at $k$-loop order, together with lower-order terms, fixes the remaining single logarithm coefficients $c_{k,k-1}$. The latter well-known RG properties constitute a crucial preliminary step of our RGOPT calculations, first requiring the precise perturbative $m$ dependence, namely the relevant coefficients including massive quarks in (4.1).

At three loops accordingly all the logarithmic coefficients $c_{3i}$, $i \le 2$ are easily determined [26] as mentioned above from lower orders and RG properties. The remaining nonlogarithmic coefficient $c_{33}$, not related to RG properties, was calculated in [43] from related three-loop quantities. In our normalization (and restricted to $N_c = 3$ for QCD) these coefficients read

$$\begin{aligned} c_{30} &= -\frac{2}{9}(81 - 2n_f), \\ c_{31} &= \frac{2}{9}(141 - 5n_f), \\ c_{32} &= \frac{1}{16}\left(52n_h + 20n_l - \frac{4406}{9} + \frac{32}{3}z_3\right), \\ c_{33} &= \frac{1}{432}\left[6185 - 768a_4 - 32\ln^2(2\ln^2(2-6z_2)) + 504z_3 + 528z_4 + (672z_3 - 750)n_h - 6n_l(32z_3 + 45)\right], \end{aligned} \quad (4.3)$$
for $n_l$ “light” (massless) and $n_h$ massive quarks, with $n_f = n_l + n_h$, $z_k = \zeta(k)$, and $a_4 = \text{Li}_4(1/2)$. In Eq. (4.1) $n_l$ and $n_h$ do not enter explicitly at one and two loops (fully described, up to unshown counterterms, by the first two graphs of Fig. 1). At three loops $n_l$ and $n_h$ enter independently only within the $L_m$ and nonlogarithmic coefficients $c_{32}$, $c_{33}$ respectively, as can be deduced from inspection of the graphs of Fig. 1. To give a more numerical savor, in particular of the $n_l$, $n_h$ dependence and relative size compared to the other (pure gauge) contributions, one has to reasonable ($10^{-6}$) accuracy:

$$ 
\begin{aligned}
c_{30} &= -18 + \frac{4}{9}n_f, \\
c_{31} &= \frac{94}{3} - \frac{10}{9}n_f, \\
c_{32} &= -29.7959 + 3.25n_h + 1.25n_l, \\
c_{33} &= 16.4566 + 0.133755n_h - 1.15925n_l.
\end{aligned} 
\quad (4.4) $$

Next at higher orders using $\beta(a_s)$ and $\gamma_m(a_s)$ to four loops [44] and $\Gamma^0(a_s)$ to five loops [32], we obtain after algebra the four-loop and five-loop exact analytical expressions of the logarithmic coefficients given in the Appendix [see Eqs. (A6) and (A8)–(A12)]. Numerically at four-loop order this reads<fnref n="4" />

$$ 
\begin{aligned}
c_{40} &= 85.5 - 5.11111n_f + 0.0740741n_f^2, \\
c_{41} &= -224.333 + 16.7037n_f - 0.246914n_f^2, \\
c_{42} &= 342.151 - 51.1008n_h - 32.1008n_l + 0.975309n_h^2 \\
&\quad + 0.308642n_l^2 + +1.28395n_hn_l, \\
c_{43} &= -375.082 + 42.6214n_h + 43.5949n_l - 0.0382074n_h^2 \\
&\quad - 0.790268n_hn_l - 0.752061n_l^2,
\end{aligned} 
\quad (4.5) $$

and at five-loop order<fnref n="5" />

$$ 
\begin{aligned}
c_{50} &= -418.95 + 42.1444n_f - 1.38519n_f^2 + 0.0148148n_f^3, \\
c_{51} &= 1469.29 - 173.995n_f + 6.01543n_f^2 - 0.0617284n_f^3, \\
c_{52} &= -3079.72 + 436.666n_f - 14.1506n_f^2 + 0.102881n_f^3 \\
&\quad + n_h(155.167 - 11.7778n_f + 0.222222n_f^2), \\
c_{53} &= 5102.45 - 852.446n_h - 843.205n_l + 61.7769n_hn_l \\
&\quad + 27.7449n_h^2 + 34.032n_l^2 - 0.344719n_h^2n_l \\
&\quad - 0.701646n_hn_l^2 + 0.00406919n_h^3 - 0.352858n_l^3, \\
c_{54} &= (n_f - 24.5)c_{44} - 617.146 + 309.613n_h + 144.324n_l \\
&\quad - 16.7381n_h^2 - 4.8565n_l^2 - 21.5946n_hn_l \\
&\quad - 0.0719093n_h^2n_l - 0.0533908n_hn_l^2 - 0.0301426n_h^3 \\
&\quad - 0.0116241n_l^3,
\end{aligned} 
\quad (4.6) $$

where we made clear the $c_{54}$ dependence upon the four-loop nonlogarithmic $c_{44}$ coefficient, not yet given explicitly at this stage as this deserves a more detailed discussion in the next subsection below. Similarly we have derived all the six-loop coefficients that are determinable exactly from RG properties: These are given in the Appendix [see Eq. (A14)].

Note that the nonlogarithmic five-loop coefficient $c_{55}$ is presently not known, and this finite contribution (before renormalization) is presumably technically very challenging to evaluate. Fortunately it does not play any role in our (five-loop) determination below since as above explained in Sec. II, only the $\ln^p[m]$, $p \ge 1$ terms contribute to the spectral density, Eq. (2.3).

## B. Exact versus approximate determinations of $c_{44}$

As just mentioned, we stress that more generally all the nonlogarithmic coefficients $c_{kk}$ in Eq. (4.1) trivially do not contribute directly to the spectral density at any order $k$. Yet these $c_{kk}$ are actually indirectly relevant, depending at which perturbative order one is performing calculations, since those coefficients enter in the next-order $c_{k+1k}$ single logarithm coefficient via RG properties, as explicitly demonstrated in Eq. (4.6) [see also Eq. (A12) in the Appendix]. The four-loop nonlogarithmic $c_{44}$ coefficient was not known until very recently; nevertheless we could derive its approximate (but dominant) contribution by exploiting other known four-loop results, as explained next. However while completing the present work, interestingly the complete $c_{44}$ has been very recently calculated [45], which allows us to perform the five-loop RGOPT analysis with a fully known $c_{54}$ coefficient.

Let us first derive our approximation for $c_{44}$ (that we will also use in the numerics below, to assess the sensitivity of our method upon such variations in the perturbative coefficients). For that purpose we exploit the relation of the condensate to another four-loop contribution as follows:

<footnote n="4">⁴We should point to a correction in $c_{43}$ here as compared with Eq. (5.13) of [26] (that was also differently normalized by an overall $4^3$ factor): This mistake, due to our previously incorrect interpretation of $n_h$ dependence from given $n_l = n_f - 1$, $n_h = 1$ results, changes $c_{43}$ by a few 0.1%, but affects our four-loop RGOPT condensate value by less than $10^{-3}$.</footnote>

<footnote n="5">⁵The authors of [32] provide the vacuum-energy anomalous dimension at five loops for both diagonal contributions of the $n_l$ massless and $n_h = 1$ massive quark, and nondiagonal contributions (i.e., quarks of different masses). It is straightforward to derive from their results the more specific case of $n_l = 0$ and $n_h(=n_f)$ degenerate quarks of mass $m$, more relevant to our calculation.</footnote>
$$ \frac{\partial}{\partial m} (\langle \bar{q}q \rangle (m)) = -\Pi_s(q^2 = 0), \qquad (4.7) $$

where $\Pi_s(q^2) \equiv i \int d^4x e^{iq \cdot x} \langle 0 | T J^s(x) J^s(0) | 0 \rangle$ is the two-point scalar correlation function (the scalar current being defined as $J^s = \bar{q}q$). This well-known relation (see e.g., [20]) is valid to all orders both at the bare and renormalized levels. The various [vector, axial, (pseudo)scalar] correlators have been investigated intensively in the literature [46,47], and up to four loops [48,49]. In particular the four-loop $\Pi_s(0)$ contribution was calculated in [49], however not incorporating the so-called singlet contributions (as those were not directly relevant to the calculation of [49]). [We recall that the singlet contributions involving two disconnected quark lines in the two-point correlators, only appear starting at three-loop order, and for $\Pi_s(0)$ at three and four loops they are nonvanishing only for massive quark contributions $\propto n_h$.] The nonsinglet four-loop non-logarithmic contribution to $\Pi_s(0)$ is given in the $\overline{\text{MS}}$ scheme in Eq. (B.1) of [49], a result that we recast here for completeness in our normalization conventions:

$$ \Pi_s^{4-\text{loop},ns}(0) = \left(\frac{3}{2\pi^2}\right) a_s^3 m^2 \left(\frac{1}{2} \bar{C}_{-1}^{(3),s} + \ln(m/\mu) \text{ terms}\right), \qquad (4.8) $$

with

$$ \begin{aligned} \bar{C}_{-1}^{(3),s} = & -325.6276432 + 16.39537650n_h + 19.76434509n_l \\ & - 1.670198265n_h^2 - 0.9856898698n_l n_h \\ & + 0.7103788267n_l^2. \end{aligned} \qquad (4.9) $$

Now at the level of the quark condensate being a one-point function, there is no distinction between “singlet” and “nonsinglet” contributions, these being all included if the condensate is calculated from basics. But if deriving the condensate using Eq. (4.7), we may explicitly separate the contributions that correspond to “singlet” or “nonsinglet” within $\Pi_s(0)$. Accordingly from a straightforward integration from Eq. (4.7) with input (4.8), we obtain the “incomplete-singlet” (IS) approximation of $c_{44}^6$:

$$ \begin{aligned} c_{44}^{JS} = & -\frac{1}{6}\bar{C}_{-1}^{(3),s} - \frac{1}{3}c_{43} \\ = & 179.2986813 - 16.93968372n_h + 0.2911021798n_h^2 \\ & - 17.82568131n_l + 0.1322905474n_l^2 \\ & + 0.4277044657n_hn_l. \end{aligned} \qquad (4.10) $$

Note that the first (dominant) term in Eq. (4.10) is the pure gauge contribution, while terms $\propto n_l, n_h$ originate from four-loop contributions with virtual massless and massive quarks.

Alternatively, the independent calculation very recently performed in [45] includes the complete contributions directly for the condensate: In the normalization of Eq. (4.1) this full $c_{44}$ reads⁷

$$ \begin{aligned} c_{44} = & 179.29868127533155 - 15.013277376448457n_h \\ & + 0.7428868214454403n_h^2 - 17.825681312474572n_l \\ & + 0.13229054734724904n_l^2 \\ & + 1.0016895879838739n_l n_h. \end{aligned} \qquad (4.11) $$

As can be seen Eq. (4.10) is fully consistent with the complete result of Eq. (4.11) (numerically within $10^{-10}$ relative accuracy) for its “nonsinglet” part (including in particular the dominant gauge contributions). Numerically the additional contributions within the full $c_{44}$ are not at all negligible at four loops: For our relevant case with no massless quarks ($n_l = 0$) and $n_h(=n_f)$ (degenerate) massive ones, Eq. (4.11) is ~4% (~7.5%) larger than (4.10), respectively for $n_f = 2$ ($n_f = 3$). In the numerics below we evidently preferably use the full Eq. (4.11), relevant for the five-loop spectral density via Eq. (4.6), but in Sec. V D we also compare results obtained with the “incomplete-singlet” approximation Eq. (4.10) in order to have a sensible estimate of the stability of five-loop RGOPT results with respect to this well-defined variation of the perturbative coefficients. We anticipate that it impacts the final condensate value roughly by a 1(2)% change of the relative magnitude of $|\langle \bar{q}q \rangle|^{1/3}$ respectively for $n_f = 2$ ($n_f = 3$).

### C. Explicitly RG-invariant condensate

One may use RG properties to define a RG-invariant renormalized condensate expression, namely that obeys the homogeneous RG Eq. (3.4), by compensating for the anomalous dimension in Eq. (4.2), as follows. The RG noninvariance of (4.1) can be perturbatively restored most simply upon considering perturbative extra finite subtraction contributions [26,36],

$$ [m\langle\bar{q}q\rangle]_{\text{inv}} \equiv m\langle\bar{q}q\rangle - S(m, a_s), \qquad (4.12) $$

where we define

$$ S(m, a_s) = \frac{3}{2\pi^2} \frac{m^4}{a_s} \sum_{k \ge 0} s_k a_s^k \qquad (4.13) $$

<footnote>$^6$Equation (4.7) implies that $c_{43}$ also enters this relation. Since $c_{43}$ is an exact contribution from the condensate, Eq. (4.10) involves both “singlet” and “nonsinglet” (from $\bar{C}_{-1}^{(3),s}$) contributions to $\Pi_s(0)$.</footnote>

<footnote>$^7$The original four-loop results of [45] combine exact analytical contributions with other (gauge) contributions known numerically but to very high accuracy of at least $10^{-74}$. Here we give for compactness the results numerically with $10^{-16}$ accuracy.</footnote>
with coefficients determined order by order by

$$ \mu \frac{d}{d\mu} S(m, a_s) = \mu \frac{d}{d\mu} m(\langle \bar{q}q \rangle) = -4m^4\Gamma^0(a_s). \quad (4.14) $$

Once having determined as above all the correct logarithmic coefficients $c_{kj}$, $j < k$ at perturbative order $k$, one may apply the first equality in Eq. (4.14), using the RG operator Eq. (3.4), to the finite expression (4.1), not separately RG invariant, to determine the subtraction function $S(m, a_s)$ uniquely. Of course, $S(m, a_s)$ actually only depends on the vacuum-energy anomalous dimension and other RG functions $\beta(a_s)$ and $\gamma_m(a_s)$, as the second equality in (4.14) shows [which is nothing but a rewriting of Eq. (4.2) above]. Note that Eq. (4.13) necessarily starts with the $s_0/a_s$ term to be consistent with RG-invariance properties. In our normalization (4.13) the exact $s_i$ expressions up to five loops are given for completeness in the Appendix [see Eq. (A16)]. Note, however, that Eq. (4.13) plays actually no role in our subsequent determination of the condensate in the present work, since $S(m, a_s)$ does not involve any $\ln(m)$ terms, so trivially it does not contribute to the spectral density. We have worked out this quantity for completeness since the expression (4.12) is nevertheless useful in other context (see e.g., Ref. [29]).

### D. Perturbative spectral density at five and six loops

From the generic perturbative expansion for the condensate, Eq. (4.1) calculating the (perturbative) spectral density formally involves calculating all logarithmic discontinuities according to Eq. (2.3). This is simply given by taking in (4.1) all nonlogarithmic terms to zero, those having obviously no discontinuities, while replacing all powers of logarithms, using $m \to i|\lambda|$ etc., as

$$ \ln^n \left( \frac{m}{\mu} \right) = \frac{1}{2^n} \ln^n \left( \frac{m^2}{\mu^2} \right) \\ \to \frac{1}{2^n 2i\pi} \left[ \ln^n \left( \frac{|\lambda|^2}{\mu^2} e^{i\pi} \right) - \ln^n \left( \frac{|\lambda|^2}{\mu^2} e^{-i\pi} \right) \right] \quad (4.15) $$

leading to the following substitution rules for the first few terms

$$ \begin{aligned} \ln\left(\frac{m}{\mu}\right) &\to \frac{1}{2}, & \ln^2\left(\frac{m}{\mu}\right) &\to L_{\lambda}, \\ \ln^3\left(\frac{m}{\mu}\right) &\to \frac{3}{2}L_{\lambda}^2 - \frac{\pi^2}{8}, & \ln^4\left(\frac{m}{\mu}\right) &\to 2L_{\lambda}^3 - \frac{\pi^2}{2}L_{\lambda}, \\ \ln^5\left(\frac{m}{\mu}\right) &\to \frac{5}{2}L_{\lambda}^4 - \frac{5\pi^2}{4}L_{\lambda}^2 + \frac{\pi^4}{32}, & & \\ \ln^6\left(\frac{m}{\mu}\right) &\to 3L_{\lambda}^5 - \frac{5\pi^2}{2}L_{\lambda}^3 + 3\frac{\pi^4}{16}L_{\lambda}, & & \end{aligned} \quad (4.16) $$

where $L_{\lambda} \equiv \ln(\lambda/\mu)$ (note the $\pi^{2k}$ terms appearing starting at order $\ln^3 m$).

We obtain in this way the perturbative spectral density up to six-loop order formally,

$$ -\rho_{\text{QCD}}^{\overline{\text{MS}}}(\lambda, a_s) = \left(\frac{3}{2\pi^2}\right) \lambda^3 \left\{ -\frac{1}{2} + 4a_s \left(L_{\lambda} - \frac{5}{12}\right) \right. \\ \left. + a_s^2 \sum_{i=1}^3 \rho_{3i} L_{\lambda}^{3-i} + a_s^3 \sum_{i=1}^4 \rho_{4i} L_{\lambda}^{4-i} \right. \\ \left. + a_s^4 \sum_{i=1}^5 \rho_{5i} L_{\lambda}^{5-i} + a_s^5 \sum_{i=1}^6 \rho_{6i} L_{\lambda}^{6-i} \right\}, \quad (4.17) $$

where the coefficients $\rho_{ki}$ for $k \ge 3$ are straightforwardly related to the $c_{ki}$ of the original condensate using (4.16) as follows:

$$ \begin{aligned} \rho_{31} &= \frac{3}{2}c_{30}, \\ \rho_{32} &= c_{31}, \\ \rho_{33} &= \frac{1}{2}c_{32} - \frac{\pi^2}{8}c_{30}, \end{aligned} \quad (4.18) $$

$$ \begin{aligned} \rho_{41} &= 2c_{40}, \\ \rho_{42} &= \frac{3}{2}c_{41}, \\ \rho_{43} &= -\frac{\pi^2}{2}c_{40} + c_{42}, \\ \rho_{44} &= -\frac{\pi^2}{8}c_{41} + \frac{1}{2}c_{43}, \end{aligned} \quad (4.19) $$

$$ \begin{aligned} \rho_{51} &= \frac{5}{2}c_{50}, \\ \rho_{52} &= 2c_{51}, \\ \rho_{53} &= -\frac{5}{4}\pi^2 c_{50} + \frac{3}{2}c_{52}, \\ \rho_{54} &= -\frac{\pi^2}{2}c_{51} + c_{53}, \\ \rho_{55} &= \frac{\pi^4}{32}c_{50} - \frac{\pi^2}{8}c_{52} + \frac{1}{2}c_{54}, \end{aligned} \quad (4.20) $$

and so on at higher (six-loop) order [see Eq. (A15) in the Appendix].

At this stage, before proceeding with RGOPT, we remark that the above (ordinary) perturbative spectral density $\rho(\lambda)$ expression for arbitrary $n_l, n_h$ in Eq. (4.17) can be useful for different purposes, independent of the RGOPT approach. For instance it should allow us to proceed at higher order the recently developed approach of Ref. [50], to fit recent lattice precise calculations of the spectral density, in order to extract $\alpha_s$.
## V. NUMERICAL RGOPT RESULTS FOR THE CONDENSATE UP TO SIX LOOPS

We are now fully equipped to proceed with the main purpose, that we recap is to find solutions of the RGOPT equations, Eqs. (3.10) and (3.8), applied to the spectral density Eq. (4.17) at successive orders, after the modifications implied by Eqs. (3.9) and (3.11). The RGOPT results up to four loops were obtained in [26] to which we refer for more details. Here we will first summarize the main steps and important features for self-containedness, before presenting in more details our new results at five and six loops. We also discuss in details the numerical impact of some controllable approximations, that will be specified, and how we accordingly estimate theoretical uncertainties of our predictions.

### A. Summary of previous results up to four loops

At one-loop order $O(1)$ for the spectral density, (3.9) only affects the first constant term $-1/2$ in Eq. (4.17): Since there is no logarithmic $L_\lambda$ contribution, one obtains for Eq. (3.10) the trivial optimized solution, $\lambda = 0$. Thus nontrivial solutions occur starting at next-to-leading (NLO) two-loop order of the modified perturbation. Accordingly at NLO the modified series reads

$$ -\rho_{\text{QCD}}^{\delta^1} = \frac{3}{2\pi^2} \lambda^3 \left( \frac{19}{58} + \frac{g}{\pi^2} \left( L_{\lambda} - \frac{5}{12} \right) \right), \quad (5.1) $$

and the OPT (3.10) and RG (3.8) equations have a unique solution given in the first lines of Tables I and II for $n_f = 2$, 3 respectively, using also (2.2). These results used the RG Eq. (3.8) at one-loop order, that gives simple analytic solutions. But since our optimized expression actually relies on exact two-loop calculations, it appears more sensible to use the RG Eq. (3.8) at the same (two-loop) order to incorporate a priori more consistently higher-order

<table><caption>TABLE I. $n_f = 2$ RGOPT results at successive orders up to four loops for the spectral parameter $\tilde{\lambda}$, $\tilde{\alpha}_S$, and RG-invariant condensate $\langle \bar{q}q \rangle_{\text{RGI}}^{1/3}$ calculated at the consistent perturbative order from (5.3). $\tilde{\Lambda}_2$ is conventionally normalized in most cases by Eq. (5.2), except in the very first line where the one-loop expression $\tilde{\Lambda} = \mu e^{-1/(2b_0 g)}$ is rather used. The corresponding scale values from Eq. (5.2) are also given in the last column.</caption><thead><tr><th>δ<sup>k</sup>, RG order</th><th>ln(λ̃<sub>μ</sub>)</th><th>α̃<sub>S</sub></th><th>(-⟨q̃q⟩<sub>RGI</sub>)<sup>1/3</sup>/λ̃<sub>2</sub></th><th>μ̃/λ̃<sub>2</sub></th></tr></thead><tbody><tr><td>δ, RG one loop</td><td>-2275/10092</td><td>87π/328 ≈ 0.83</td><td>0.996</td><td>2.2</td></tr><tr><td>δ, RG two loop</td><td>-0.45</td><td>0.480</td><td>0.821</td><td>2.8</td></tr><tr><td>δ<sup>2</sup>, RG two loop</td><td>-0.686</td><td>0.483</td><td>0.792</td><td>2.797</td></tr><tr><td>δ<sup>2</sup>, RG three loop</td><td>-0.703</td><td>0.430</td><td>0.783</td><td>3.104</td></tr><tr><td>δ<sup>3</sup>, RG three loop</td><td>-0.83895</td><td>0.40522</td><td>0.77428</td><td>3.308</td></tr><tr><td>δ<sup>3</sup>, RG four loop</td><td>-0.82164</td><td>0.39071</td><td>0.77247</td><td>3.448</td></tr></tbody></table>

<table><caption>TABLE II. Same caption as Table I for $n_f = 3$.</caption><thead><tr><th>δ<sup>k</sup> order</th><th>ln(λ̃<sub>μ</sub>)</th><th>α̃<sub>S</sub></th><th>(-⟨q̃q⟩<sub>RGI</sub>)<sup>1/3</sup>/λ̃<sub>3</sub></th><th>μ̃/λ̃<sub>3</sub></th></tr></thead><tbody><tr><td>δ, RG one loop</td><td>-283/972</td><td>27π/104 ≈ 0.82</td><td>0.987</td><td>2.35</td></tr><tr><td>δ, RG two loop</td><td>-0.56</td><td>0.474</td><td>0.789</td><td>3.06</td></tr><tr><td>δ<sup>2</sup>, RG two loop</td><td>-0.766</td><td>0.493</td><td>0.772</td><td>2.942</td></tr><tr><td>δ<sup>2</sup>, RG three loop</td><td>-0.788</td><td>0.444</td><td>0.766</td><td>3.273</td></tr><tr><td>δ<sup>3</sup>, RG three loop</td><td>-0.97402</td><td>0.41367</td><td>0.74377</td><td>3.547</td></tr><tr><td>δ<sup>3</sup>, RG four loop</td><td>-0.96506</td><td>0.39906</td><td>0.74232</td><td>3.709</td></tr></tbody></table>

effects. Doing this gives the results in the second lines of Tables I and II for $n_f = 2, 3$. Those results, to be considered more accurate, show a substantial decrease of the optimal coupling $\alpha_S$ to a more perturbative value with respect to the results using the one-loop RG equation.

At higher orders the precise numbers obtained for the condensate also depend on the specific definition of the $\tilde{\Lambda}$ reference scale, which is generally perturbative and a matter of convention to some extent. The numbers in the first lines of Table I were obtained using the simpler one-loop form, $\tilde{\Lambda} = \mu e^{-1/(2b_0 g)}$, consistent with the one-loop RG equation used. Next, when comparing below with other determinations of the condensate, we use conventionally a four-loop definition of $\tilde{\Lambda}$ (see, e.g., [51]), in agreement with most other past determination conventions. Except, at five-loop order, we obviously adopt the more consistent five-loop perturbative definition of $\tilde{\Lambda}$. The four-loop QCD scale $\tilde{\Lambda}$ expression reads in our normalizations [where $g = 4\pi\alpha_S(\mu)$],

$$ \tilde{\Lambda}_{n_f}^{4-\text{loop}}(g) \equiv \mu e^{-\frac{1}{2b_0 g}} (b_0 g)^{-\frac{b_1}{2b_0^2}} \exp \left[ -\frac{g}{2b_0} \cdot \left( \left( \frac{b_2}{b_0} - \frac{b_1^2}{b_0^2} \right) + \left( \frac{b_1^3}{2b_0^3} - \frac{b_1 b_2}{b_0^2} + \frac{b_3}{2b_0} \right) g \right) \right], \quad (5.2) $$

with a straightforward generalization upon including the five-loop coefficient $b_4$. In Tables I and II we actually give for convenience the value of the scale-invariant condensate $\langle \bar{q}q \rangle_{\text{RGI}}$, which can be more appropriately compared between different perturbative orders. It is defined in our normalization as

$$ \begin{aligned} \langle \bar{q}q \rangle_{\text{RGI}} &\equiv \langle \bar{q}q \rangle(\mu) \exp \left[ \int dg \frac{\gamma_m(g)}{\beta(g)} \right] \\ &= \langle \bar{q}q \rangle(\mu) (2b_0 g)^{\frac{\gamma_0}{2b_0}} \\ &\quad \times \left( 1 + \left( \frac{\gamma_1}{2b_0} - \frac{\gamma_0 b_1}{2b_0^2} \right) g + \mathcal{O}(g^2) \right) \end{aligned} \quad (5.3) $$

where higher-order terms not shown here are easily derived, since only depending on the RG coefficients $b_i$, $\gamma_i$ known up to five loops. We remark, however, that our RGOPT optimization also fixes a scale simply obtained from using
<figure><img src="image_4.png" /><figcaption>FIG. 2. The different branch solutions $L_\lambda^{\text{RG}}(\alpha_s, k)$ (solid curves) and $L_\lambda^{\text{OPT}}(\alpha_s, k)$ (dashed curves), $k = 1, 2, 3$, at $\delta^3$ (four-loop) order. Left: $n_f = 3$. Right: $n_f = 2$. Distinct branches appear to join within some $\alpha_s$ range where they actually become complex conjugated, since only their real parts are plotted within this range.</figcaption></figure>

Eq. (5.2) (or its lower-order equivalent) for $\bar{\Lambda}(\tilde{g})$, that is given indicatively in Tables I and II. We stress that the optimal coupling $\tilde{\alpha}_S$ and corresponding optimal scale $\tilde{\mu}$, or the optimal spectral parameter $\tilde{\lambda}$, are to be considered intermediate values with no universal physical interpretation, since their precise obtained values depend on the physical quantity being optimized. The physically meaningful result is obtained when inserting $\tilde{\alpha}_S$ and $\tilde{\lambda}$ within the quantity being optimized, here $\rho(\lambda, \alpha_S)$. (This feature is quite general in optimization procedures: The values of the optimization parameters for a given physical quantity should not in general be used to evaluate another physical quantity.)

At three-loop, $a_s^2$, order, the $n_f$ dependence appears explicitly within the perturbative coefficients of the spectral density; see Fig. 1 and the last $a_s^2$ coefficient in Eq. (4.17). There occur two real solutions for $\tilde{L}_\lambda$, $\tilde{\alpha}_S$, but the selection of the unique physical solution is unambiguous since only one is clearly compatible with AF behavior for $g \to 0$, $\ln(\tilde{\lambda}/\mu) \simeq -d_k/(2b_0g) + \mathcal{O}(1)$ with $d_k = \mathcal{O}(1)$, both for the RG and OPT equations. In contrast the other real solution has for $g \to 0$ a coefficient of opposite sign to AF, and gives $\ln \tilde{\lambda}/\mu > 0$, which is incompatible with perturbativity, since we expect $\mu \gg \tilde{\lambda}$ similar to the perturbative range $\mu \gg \tilde{m} \sim \bar{\Lambda}$ for the original expansion with mass dependence. As stressed above in Sec. III the occurrence of an essentially unique solution with the correct AF-matching behavior at successive orders is a crucial feature of RGOPT, as will be illustrated further below.

At three and four loops the RGOPT results for $n_f = 2, 3$ are specified in Tables I and II respectively.<fnref n="8" /> As indicated in each case we compare results obtained when using first the

RG Eq. (3.8) truncated at lower order, and next taking the full RG equation at the same three- or four-loop order respectively, incorporating more higher-order dependence. At four loops the raw optimization results actually give several real solutions for $\tilde{\lambda}, \tilde{\alpha}_S$ but there are no possible ambiguities since once more all solutions are eliminated from the AF-matching requirement, except a single one, with $\tilde{\alpha}_S > 0$ and $\tilde{L}_\lambda < 0$ as expected.

One observes a further decrease of the optimal coupling $\tilde{\alpha}_S$ from three to four loops to more perturbative values, as well as the corresponding decrease of $\tilde{L}_\lambda$, meaning that $\tilde{\mu}$ is also larger. The stabilization/convergence of the results is clear for the scale-invariant condensate $\langle \bar{q}q \rangle_{\text{RGI}}$ given in Tables I and II, which at four-loop order has almost no variation upon RG equation truncations.<fnref n="9" /> Note that the optimal values $\tilde{\alpha}_S$ decrease substantially with increasing orders as compared to the lowest nontrivial order result above, thus indicating more perturbatively reliable results; moreover $\tilde{\alpha}_S$ appears to somehow stabilize at three and four loops. Notice also that compared with the more than 10% change in $\tilde{\alpha}_S$ upon going from two to three loops, the final physical condensate value only varies by 0.25%, showing a strong stability. Also, while $\langle \bar{q}q \rangle^{1/3}/\bar{\Lambda}$ changes by about 20% compared to the crude two-loop result in the first lines of Tables I and II it stabilizes rapidly at higher orders showing a *posteriori* that the first nontrivial two-loop result seems already a quite realistic value. This stability at only NLO is a welcome feature for the usefulness of the RGOPT. A similar behavior was observed when optimizing the pion decay constant in [25].

In Fig. 2 we illustrate for $n_f = 3$ (resp. left) and $n_f = 2$ (resp. right) the different RG and OPT branches obtained at

<footnote n="8">⁸Since all the perturbative coefficients are known exactly at four loops, or to very high accuracy at five loops, our optimized results at a given order are in principle obtained to high accuracy. But in Tables I and II (and similarly at higher orders below) we give results to an accuracy largely sufficient for our purpose, given the extra uncertainties that will be discussed below.</footnote>

<footnote n="9">⁹In our numerical analysis below we use for convenience the exponentiated form of the RG-invariant factor as in (5.3), but note that the relative difference with the fully perturbatively expanded one is less than $10^{-3}$ for all considered optimized coupling values.</footnote>
four loops, respectively as $L_{\lambda}^{\text{RG}}(\alpha_s, k)$, $L_{\lambda}^{\text{OPT}}(\alpha_s, k)$, where the number of solutions are at most $k = 1..3$ at four loops. The clearest situation is the one for $n_f = 3$, where the two AF-matching RG and OPT branches (namely the two curves with the lowest $L_{\lambda}$ values for any $\alpha_s > 0$ in Fig. 2) are real for any $\alpha_s > 0$ and intersect at a unique value, that determines unambiguously the solution; compare the last line in Table II. Similar properties hold for all considered cases at lower orders. Next in Fig. 2 (right) for $n_f = 2$, two of the OPT branches inconveniently become complex (conjugate) valued within the range $0.25 \lesssim \alpha_s \lesssim 0.33$ (so that their real parts shown here appear joined). Nevertheless one can still unambiguously select the correct AF-matching OPT branch, that is the one intersecting with the AF-matching RG branch, so that the correct solution is again unique.

## B. Five-loop and six-loop results

Up to four loops, all the perturbative coefficients and RG quantities entering our evaluation have been known exactly for some time. Thanks to the recently calculated five-loop vacuum anomalous dimension [32], Eq. (A5), and the very recent complete calculation [45] of the four-loop nonlogarithmic coefficient $c_{44}(n_f)$ given in Eq. (4.11), all the relevant perturbative coefficients needed at five loops are exactly available for the spectral density, Eq. (4.17). Thus we can extend our evaluation for the physically relevant $n_f = 2, 3$ values to five-loop order, correspondingly including up to five-loop contributions in the RG $\beta$ and $\gamma_m$ functions within the optimization, after performing consistently the $\delta$ expansion in Eq. (3.9) to order $\delta^4$. As mentioned above it is also useful to estimate the sensitivity of our results to the well-defined approximation Eq. (4.10) neglecting in $c_{44}$ the (subdominant) four-loop singlet contributions. Furthermore, higher-order coefficients are (partly) determinable solely from perturbative RG invariance, a feature that we can exploit to consider also (approximate) six-loop results (RGOPT order $\delta^5$). More precisely all the presently known five-loop RG coefficients, together with the complete four-loop coefficients, allow us to determine exactly the six-loop coefficients $c_{6k}$, $0 \le k \le 4$ of $\ln^{6-k}(m/\mu)$ of Eq. (4.1). While the single logarithmic term ($k=5$) would need the presently unknown six-loop vacuum-energy anomalous dimension as well as the five-loop nonlogarithmic coefficient $c_{55}$. The explicit expressions of the $c_{6k}$ are given in Eq. (A14) in the Appendix. Consequently from Eq. (4.16) all six-loop logarithmic terms, $\rho_{6k} L_{\lambda}^{6-k}$, $k = 1, 5$ of the spectral density in Eq. (4.17) are exactly predicted [see Eq. (A15)], except for its last unknown nonlogarithmic coefficient. As we will examine below, the six-loop results, although being approximate, are quite important to assess a more reliable determination of the condensate, due to the occurrence of rather unwelcome instabilities for the strictly five-loop results.

## 1. Five-loop and six-loop $n_f = 2$ results

We examine now in some details our procedure and results for $n_f = 2$, with quite similar features given more briefly below for $n_f = 3$, except when important differences need to be mentioned. For $n_f = 2$, at five loops we obtain one real solution that appears at first sight the closest to the lower (four-loop) results, namely, $L_{\tilde{\lambda}} = -0.5699$, $\tilde{\alpha}_s = 0.5963$, which gives $\langle \bar{q}q \rangle_{\text{RGI}}^{1/3} \approx -0.863 \bar{\Lambda}_2$ obtained using the RG equation at five loops. (Very close results are obtained if using instead the RG equation at four loops.) Without further inquiries one would conclude from this result that the five-loop RGOPT produces an anomalously large shift of the condensate value, as compared with the seemingly well-stabilized three- and four-loop results $\sim - (0.78 - 0.77)$ in Table I. A directly related issue is the anomalously large optimized coupling that corresponds to this solution, $\alpha_s \sim 0.6$, in contrast with the regularly decreasing coupling obtained at increasing orders up to four loops, in Tables I and II.

However, upon applying our general criteria to select the correct solutions, a more careful examination shows that this solution cannot be correct, since it is not sitting on the perturbative AF-matching branch, in contrast to what occurs systematically at lower orders. This feature can be checked rather easily by perturbatively expanding at first order the four different branch solutions for RG and OPT equations respectively, that both give quartic equations in $L_{\lambda} \equiv \ln(\lambda/\mu)$ at five loops [thus respectively giving $L_{\tilde{\lambda}}^{\text{RG}}(\alpha_s, k)$, $L_{\tilde{\lambda}}^{\text{OPT}}(\alpha_s, k)$ with $k = 1, \dots, 4$], and examining which one(s) exhibit the perturbative AF-matching behavior, and whether the latter are matching the optimized $L_{\tilde{\lambda}}$, $\tilde{\alpha}_s$ values obtained at the intersecting solution(s). Equivalently it can be seen more pictorially in Fig. 3 illustrating the different branches and (some of) their intersecting solutions (these branches are shown in a somewhat restricted but physically relevant range of $L_{\lambda}$ and $\alpha_s$): In contrast with the four-loop results in Fig. 2, the RG branches now also become complex, similar to the OPT branches, within a rather important $\alpha_s$ range $0.27 \lesssim \alpha_s \lesssim 0.56$, and the only real intersection occurring at $\alpha_s \approx 0.596$, $L_{\tilde{\lambda}} \approx -0.57$ (visible near the top right of Fig. 3) sits on RG and OPT branches that are not linked to the correct AF behavior.

Therefore, at five-loop order the RG and OPT AF-matching branch do not have a real-valued common intersection: a feature which somewhat complicates our investigation as compared with lower orders. This large perturbation sufficiently destabilizing the regular trend observed at lower orders to suppress real AF-matching solutions, has two distinct, clearly identified origins. The first feature (but having a rather moderate impact on final results) is that the five-loop vacuum anomalous dimension coefficient [32] is much larger relative to lower orders, moreover varying very much with $n_f$ [compare $\Gamma_4^0$ in
<figure><img src="image_5.png" /><figcaption>FIG. 3. Some of the different RG (solid curves) and OPT (dashed curves) branches, respectively $L_{\lambda}^{\text{RG}}(\alpha_s)$, $L_{\lambda}^{\text{OPT}}(\alpha_s)$ at $\delta^4$ (five-loop) order for $n_f = 2$. Distinct branches appear to join within some $\alpha_s$ range where they actually become complex conjugate, since only their real parts are plotted within this range.</figcaption></figure>

Eq. (A5) with $\Gamma_3^0$ in Eq. (A4) in the Appendix]. [In contrast, as illustrated below, going from the four-loop to the five-loop $\beta$ function for the RG equation has a very modest impact, which can be traced to the moderate numerical changes of adding five-loop RG coefficients to $\beta(\alpha_s)$, $\gamma_m(\alpha_s)$ functions.] Note that $\Gamma_4^0$ enters the five-loop condensate $L_m$ coefficient $c_{54}(n_f)$, but it is not the sole contribution: The net effect from $\Gamma_4^0$ is typically that $|c_{54}(n_f)/c_{50}(n_f)| \sim 10$ roughly, which may be compared qualitatively with the similar four-loop quantities giving $|c_{43}/c_{40}| \sim 3$. But the second feature, that upon inspection happens to be the principal reason why the AF-matching solution is pushed into the complex domain, is that the discontinuities from Eq. (4.15) entail also a relatively large term $\sim \pi^4$ appearing for the first time at five-loop order in $\rho(\lambda, \alpha_s)$, within the nonlogarithmic coefficient. More precisely, it is the first term of Eq. (4.20) modifying the relevant original perturbative coefficient, $c_{54}(n_f)$ in Eq. (4.1) by about 60%, while all other coefficients are more moderately affected by the discontinuity contributions. Simply ignoring this contribution would be clearly inconsistent, and we will examine below how to better circumvent those problems.

Note in Fig. 3 that the AF-matching RG and OPT branches have not disappeared but just became complex (conjugate) valued within a certain $\alpha_s$ range, rather unfortunately where the sought intersecting solution is expected. Indeed one can determine precisely the complex-conjugated solution that sits on the AF-matching branch: Using the four-loop RG equation, we obtain $L_{\tilde{\lambda}} \simeq -0.778 \pm 0.303i$, $\tilde{\alpha}_s \simeq 0.358 \pm 0.0537i$, that gives for the (RG-invariant) condensate $\simeq -(0.801 \pm 0.0195i)\tilde{\Lambda}(2)$ (see Table III for more details). Accordingly the correct AF-matching branch, although complex valued in the relevant range, happens to give a corresponding condensate value with a small imaginary part, and with a real part in smoother continuity with the

four-loop real solution. Also the corresponding (real part of the) optimal coupling $\tilde{\alpha}_s$ is more reasonably smaller than for the (wrong) naive real solution above. Very similar results are obtained if using rather the five-loop RG equation (see Table III).

At this stage without further investigation one may just take the real part as the physically relevant result, and interpret the imaginary parts as a rough estimate of the theoretical uncertainties of the results (although this is presumably not the best possible prescription to estimate the intrinsic uncertainties). But given that the unwelcome occurrence of nonreal solutions is only a consequence of solving *exactly* the RG and OPT polynomial equations in $L_{\lambda}$, and that it is seemingly not far from a real solution at five loops, one can more appropriately attempt to recover real solutions by a variant of the procedure. Accordingly a first possibility is simply to (perturbatively) approximate the sought optimized solutions at five loops. Alternatively another possibility is to proceed to next (six-loop) order: At least this is possible in the approximation of neglecting the nonlogarithmic six-loop coefficient being the only contribution not presently derivable from already known lower-order results (as explained above at the beginning of Sec. V B). Let us examine in turn those two possibilities.

## 2. Perturbatively truncated five-loop RG solutions

At five loops, instead of solving exactly the relevant RG and/or OPT optimization Eqs. (3.8) and (3.10), one can consider more perturbative approximations, as long as those remain consistent with the original perturbative order considered. Indeed the RG Eq. (3.8) generates terms of formally higher order than five loops: More precisely it is easy to see that at five loops Eq. (3.8) acting on the five-loop ($\alpha_s^4$) spectral density Eq. (4.17) involves up to $\alpha_s^9$ terms, due to the highest five-loop RG contributions $\propto b_4\alpha_s^6$, $\gamma_4\alpha_s^5$ respectively. But $b_4$, $\gamma_4$ appear first at order $\alpha_s^5$, $\alpha_s^4$ respectively. Accordingly a presumably sensible procedure is to truncate [24] the RG equation suppressing higher-order terms in $\alpha_s$ until possibly recovering a real common RG and OPT solution. At the same time if suppressing too many higher-order terms one loses the consistency with the RG content required at a given (here four- or five-loop) order. A similar reasoning shows that the next-order six-loop RG coefficients, $b_5$, $\gamma_5$ (presently not known), would enter first respectively the $\alpha_s^6$, $\alpha_s^5$ coefficients of the RG equation.

Therefore it appears sensible to truncate any $\alpha_s^k$, $k \ge 6$ in the result of Eq. (3.8), that would be anyway affected by presently unknown higher orders. Further truncating the $\alpha_s^5$ term implies, however, losing any dependence from the five-loop $b_4$ (while it still involves the five-loop $\gamma_4$ one). Accordingly we found it instructive to consider the effects of successive truncations, progressively suppressing the highest $\alpha_s^9$ down to $\alpha_s^6$ (or even possibly $\alpha_s^5$) terms and
comparing. This is done below, with all results compiled in Table III obtained by optimizing the spectral density $\rho(\lambda, \alpha_S)$ and keeping only the AF-matching branch solution (unique at a given order).

From the $n_f = 2$ results of Table III, at five loops it appears not that easy to recover real solutions: Upon truncating terms progressively starting from highest-order ones, the results do not change much at first, although there is a slow but clear decrease of the corresponding imaginary parts. Also, despite the not small $\text{Im}[L_\lambda]$ values, the resulting condensate has much smaller imaginary parts, and real parts remain very stable, differing relatively only by $O(10^{-3})$ for the different truncations. Similarly, for all cases there are tiny differences between the results using the four-loop or five-loop RG equation. The $\tilde{\alpha}_S$ value is also more reasonably perturbative, and close to the real four-loop results of Table I. Truncating maximally the RG equation (namely by all terms $\alpha_S^{k \ge 6}$, but that still involves all the five-loop RG coefficients), the correct (AF-matching) solution has a tiny imaginary part, so that its real part may be considered reliable, giving $\langle \bar{q}q \rangle_{\text{RGI}}^{1/3} \approx -(0.800 \pm 0.0005i)\tilde{\Lambda}_2$. Note that if further truncating the RG equation, one not only loses the consistent RG content at five-loop order, but the corresponding RG equation no longer gives any AF-matching branch.

### 3. Approximate (partial) six-loop RGOPT

As sketched above, the second alternative is to proceed at next (six-loop) order of Eq. (4.17), with the $\ln 6^{-k}(\lambda)$ coefficients given explicitly in the Appendix [see Eqs. (A14) and (A15)]. The motivation, apart simply from the fact that most of the six-loop coefficients are readily exploitable from RG properties, is that the discontinuities (4.15) entail additional contributions $\propto \pi^4$ [see the last terms in Eq. (4.16)], that tend to partially balance the instability triggered by $\pi^4$ discontinuity terms appearing first at five-loop order. At this point it is worth remarking that such features are generically expected from Eq. (4.15): Typically in [26] we have calculated the spectral density for the Gross-Neveu (GN) $O(N)$ model [52] in the large-$N$ limit, to very high perturbative orders, that exhibits a clear pattern. The RGOPT solutions at increasing orders converge slowly toward the exact result (known for the GN model), those solutions being destabilized each time novel $\pi^{2k}$ contributions appear first, at increasing orders.<fnref n="10" /> However, if keeping only fixed $\pi^{2k}$ terms (namely discarding $\pi^{2k+2}$, etc., terms appearing at higher orders),

remarkably at sufficiently high fixed order all the $\pi^{2k}$ terms cancel, and the exact GN spectral density is obtained [26].

For the QCD spectral density such exact cancellations are not expected; moreover obviously we are quite limited in trying to reach still higher orders. But inspired from these properties it is worth comparing two available successive orders (five and six loops), that actually rely on the same five-loop RG content, since as we recall, five-loop RG properties predict most of the six-loop coefficients of $\rho(\lambda)$ [all except the nonlogarithmic one, $\rho_{66}$ in Eq. (4.17)]. Accordingly one should keep in mind that it remains an approximation to the complete six-loop results, since $\rho_{66}$ involves the presently unknown six-loop vacuum anomalous dimension and the five-loop nonlogarithmic coefficient $c_{55}$. Therefore we simply set $\rho_{66}$ to zero in our numerics [neglecting also consistently the other nonlogarithmic contributions generated at six loops from the discontinuities (4.16)] We will argue below that this approximation should moderately deviate from the complete six-loop results.

For $n_f = 2$ the corresponding partial six-loop RGOPT results are given in the last two lines of Table III, also considering the (maximal) RG-consistent truncation. As one can see a real solution is recovered at six loops; moreover the two AF-matching RG and OPT branches remain real for all the physically relevant $\alpha_S$ range, and their intersection occurs for a substantially smaller $\tilde{\alpha}_S$ value as compared to five loops. This is illustrated also in Fig. 4, zooming on the RG and OPT branches in the relevant range of $L_\lambda, \alpha_S$, which looks qualitatively more similar to the four-loop $n_f = 3$ case. It is striking that the resulting condensate value is much closer to the four-loop results; that is not a numerical accident but is more essentially the effect of partially balancing at six loops the instability from the large $\pi^4$ terms occurring first at five loops.

### 4. Summary of $n_f = 2$ results

As a tentative summary of the previous $n_f = 2$ investigation:

(i) At five loops, the impact of both large five-loop vacuum-energy anomalous dimensions and (more importantly) the first occurrence of $\pi^4$ terms from (4.15) are strong enough to destabilize the regular features observed at lower orders up to four loops. Consequently one fails to obtain a strictly real AF-matching solution. Yet the five-loop results from successive truncations of nonmandatory higher-order terms in the RG equation are very consistent, reflecting a good stability. Also the imaginary parts are small enough (especially for the maximal truncation of $\alpha_S^{k \ge 6}$; see Table III) and can be included within the theoretical uncertainties.

(ii) Next, going to six loops restores a real unique AF-matching solution, that results from a partial balance

<footnote n="10">¹⁰The GN spectral density exhibits at low orders even a more pronounced destabilization than for QCD, because the (large-$N$) basic perturbative expansion of $\langle \bar{q}q \rangle_{\text{GN}}(m)$ in the $\overline{\text{MS}}$ scheme has vanishing nonlogarithmic coefficients. Therefore, relative to zero, the large contributions generated by (4.15) within $\rho(\lambda)$ are maximally destabilizing corrections.</footnote>
<table><caption>TABLE III. $n_f = 2$ results at RGOPT $\delta^4$ (five loops) and partial $\delta^5$ (six loops) for the (RG-invariant) condensate $\langle \bar{q}q \rangle_{\text{RGI}}^{1/3}/\bar{\Lambda}_2$ and the corresponding optimal values of the spectral parameter $\tilde{\lambda}$, coupling $\tilde{\alpha}_S$, and scale $\tilde{\mu}$, from optimizing $\rho(\lambda, \alpha_S)$ with Eqs. (3.10) and (3.8). We compare the results of (perturbatively consistent) successive RG equation truncations. $\bar{\Lambda}_2$ is normalized by Eq. (5.2) when the four-loop RG equation is used or by its five-loop extension when the five-loop RG equation is used.</caption><thead><tr><th>RGOPT [$\rho(\lambda, \alpha_S)$] $\delta^k$, RG order</th><th>ln<sub>μ</sub><sup>λ̃</sup></th><th>α̃<sub>S</sub></th><th>-(q̃q)<sub>RGI</sub><sup>1/3</sup>/λ̃<sub>2</sub></th><th>μ̃/λ̃<sub>2</sub></th></tr></thead><tbody><tr><td>δ<sup>4</sup>, RG four loop (full)</td><td>-0.77785 ± 0.30316i</td><td>0.35785 ± 0.053706i</td><td>0.80084 ± 0.019516i</td><td>3.6418</td></tr><tr><td>δ<sup>4</sup>, RG five loop (full)</td><td>-0.78815 ± 0.31499i</td><td>0.35181 ± 0.0501i</td><td>0.80104 ± 0.019553i</td><td>3.7403</td></tr><tr><td>δ<sup>4</sup>, RG five loop (α<sub>S</sub><sup>k≥9</sup>)</td><td>-0.80175 ± 0.3081i</td><td>0.35384 ± 0.043583i</td><td>0.80101 ± 0.019379i</td><td>3.7596</td></tr><tr><td>δ<sup>4</sup>, RG four loop (α<sub>S</sub><sup>k≥8</sup>)</td><td>-0.81664 ± 0.29248i</td><td>0.35958 ± 0.036002i</td><td>0.80033 ± 0.018530i</td><td>3.7292</td></tr><tr><td>δ<sup>4</sup>, RG five loop (α<sub>S</sub><sup>k≥8</sup>)</td><td>-0.8151 ± 0.28951i</td><td>0.36097 ± 0.036464i</td><td>0.8009 ± 0.018641i</td><td>3.7094</td></tr><tr><td>δ<sup>4</sup>, RG four loop (α<sub>S</sub><sup>k≥7</sup>)</td><td>-0.83759 ± 0.25140i</td><td>0.37539 ± 0.024696i</td><td>0.80048 ± 0.015453i</td><td>3.5816</td></tr><tr><td>δ<sup>4</sup>, RG five loop (α<sub>S</sub><sup>k≥7</sup>)</td><td>-0.83618 ± 0.250i</td><td>0.37603 ± 0.025177i</td><td>0.80121 ± 0.01559i</td><td>3.5728</td></tr><tr><td>δ<sup>4</sup>, RG four loop (α<sub>S</sub><sup>k≥6</sup>)</td><td>-0.93852 ± 0.17354i</td><td>0.40216 ± 0.01192i</td><td>0.79914 ± 0.0004024i</td><td>3.3298</td></tr><tr><td>δ<sup>4</sup>, RG five loop (α<sub>S</sub><sup>k≥6</sup>)</td><td>-0.93840 ± 0.17332i</td><td>0.40223 ± 0.011871i</td><td>0.79994 ± 0.0004919i</td><td>3.3292</td></tr><tr><td>δ<sup>5</sup>, RG five loop (full)</td><td>-1.0846</td><td>0.32689</td><td>0.77133</td><td>4.3737</td></tr><tr><td>δ<sup>5</sup>, RG five loop (α<sub>S</sub><sup>k≥6</sup>)</td><td>-1.1422</td><td>0.33778</td><td>0.77260</td><td>4.1671</td></tr></tbody></table>

of the destabilizing $\pi^4$ terms. This solution has very regular properties and happens to be very close to the four-loop results.

These properties are more generically confirmed from comparison with the other relevant values $n_f = 3$, or $n_f = 0$, as illustrated next.

### C. Five- and six-loop $n_f = 3$ results

For $n_f = 3$ at five loops, very similar to $n_f = 2$ there is one real RGOPT solution appearing at first the closest to the four-loop real results, using the RG equation at five loops: $L_{\tilde{\lambda}} \approx -0.693411$, $\tilde{\alpha}_S \approx 0.598$. This gives for the RG-invariant condensate, $\approx -0.813\bar{\Lambda}_3$; thus a large shift from the four-loop results of Table II. But again upon examining the AF branches these do not match this real solution, and

<figure><img src="image_6.png" /><figcaption>FIG. 4. Some of the relevant RG (solid curves) and OPT (dashed curves) branches, respectively $L_{\tilde{\lambda}}^{\text{RG}}(\alpha_S)$, $L_{\tilde{\lambda}}^{\text{OPT}}(\alpha_S)$ at $\delta^5$ (six-loop) order for $n_f = 2$.</figcaption></figure>

the correct AF-matching but complex-valued branch gives a more reasonable result, with small imaginary parts and a real part closer to the four-loop real results (see Table IV).

Similar to $n_f = 2$ we have performed a systematic analysis of all possible RG-consistent truncations. We illustrate in Table IV the more relevant results, omitting intermediate case details. Overall the behavior is quite similar to $n_f = 2$: At five loops one fails to recover strictly real AF-matching solutions, but the maximal truncation (discarding $\alpha_S^{k\ge 6}$, not losing the five-loop RG content) gives very small imaginary parts, and we will take the real part and add appropriate uncertainties in our final estimate.

Next, similar to the $n_f = 2$ results, at six-loop order one recovers a real AF-matching solution given in the last two lines in Table IV, with very regular properties and close to the four-loop results of Table II.

### D. Impact of approximated five-loop contributions

We now consider the approximation defined in Sec. IV B and relevant for $n_f = 2, 3$, of using for the four-loop non-logarithmic coefficient $c_{44}$ our expression in Eq. (4.10) derived from the related nonsinglet four-loop scalar two-point correlator [49]. We recall that at the level of the optimized spectral density, this affects results only via the five-loop single logarithmic coefficient $c_{54}$. Since the previous results in Sec. V B including the very recently determined [45] exact $c_{44}$ coefficient Eq. (4.11) are accordingly more complete, we will not include the variations resulting from this approximation within our uncertainty estimates. Nevertheless, given that the more exact five-loop results above are somewhat prevented by
<table><caption>TABLE IV. $n_f = 3$ results at five and six loops; same caption as Table III.</caption><thead><tr><th>RGOPT [$\rho(\lambda, \alpha_S)$] $\delta^k$, RG order</th><th>ln<sup>2</sup><sub>μ</sub></th><th>α<sub>S</sub></th><th>-⟨q̃q⟩<sup>1/3</sup><sub>RGI</sub> / λ<sub>3</sub></th><th>μ / λ<sub>3</sub></th></tr></thead><tbody><tr><td>δ<sup>4</sup>, RG four loop (full)</td><td>-0.92148 ± 0.29624i</td><td>0.36925 ± 0.051547i</td><td>0.76927 ± 0.016561i</td><td>3.9124</td></tr><tr><td>δ<sup>4</sup>, RG five loop (full)</td><td>-0.93013 ± 0.3038i</td><td>0.36465 ± 0.04859i</td><td>0.76838 ± 0.016134i</td><td>3.9911</td></tr><tr><td>δ<sup>4</sup>, RG four loop (α<sub>S</sub><sup>1/2</sup>)</td><td>-1.064 ± 0.13845i</td><td>0.42391 ± 0.012618i</td><td>0.77089 ± 0.0024578i</td><td>3.4384</td></tr><tr><td>δ<sup>4</sup>, RG five loop (α<sub>S</sub><sup>1/2</sup>)</td><td>-1.0639 ± 0.13833i</td><td>0.42394 ± 0.012593i</td><td>0.76975 ± 0.0023702i</td><td>3.4297</td></tr><tr><td>δ<sup>5</sup>, RG five loop (full)</td><td>-1.2340</td><td>0.33863</td><td>0.74042</td><td>4.6923</td></tr><tr><td>δ<sup>5</sup>, RG five loop (α<sub>S</sub><sup>1/2</sup>)</td><td>-1.2618</td><td>0.3449</td><td>0.74076</td><td>4.5578</td></tr></tbody></table>

instabilities from producing real solutions, it is instructive to study their sensitivity upon such a well-defined approximation. The corresponding results are shown in Table V for $n_f = 2$ and $n_f = 3$ at five and six loops.

Very similar to the results obtained with the full $c_{44}$, at five loops the RGOPT gives nonreal AF-matching solutions, but with small imaginary parts. Accordingly in Table V the results are not too different from the ones using the exact $c_{44}$ in Tables III and IV, except that the imaginary parts are somewhat smaller. For $n_f = 3$, real AF-matching solutions are recovered upon maximal truncations consistent with RG at five loops. At six loops real AF-matching solutions are also recovered, and these differ from the ones with exact $c_{44}$ in Tables III and IV by ~1% (~2%) lower in $\bar{\Lambda}$ units for $n_f = 2$ ($n_f = 3$) respectively. All these features are consistent with the fact that the loss of a real solution at five loops is essentially due to the occurrence of relatively large $\pi^4$ terms, while the ~4% (~7%) for $n_f = 2$

($n_f = 3$) decrease in the approximated $c_{44}$ Eq. (4.10), as compared with the complete one Eq. (4.11), has a more moderate impact.

We conclude that the approximation neglecting the four-loop singlet contributions within $c_{44}$ produces a change in the final condensate magnitude $|\langle\bar{q}q\rangle|^{1/3}$ that is ~1% (2%) smaller in magnitude respectively for $n_f = 2$ ($n_f = 3$), which again reflects a good overall stability.

### E. The condensate in the quenched approximation

One can also easily extend our calculations formally for $n_f = 0$: Actually a quark of mass $m$ setting the overall scale in Eq. (4.1) “dressed” at higher orders by pure gauge interactions is still understood in this case, and the perturbative removal of the quarks entering at three loops in Fig. 1 and higher orders may be viewed as a perturbative analog of the “quenched” approximation, which has its own theoretical interest, and can be compared with lattice

<table><caption>TABLE V. $n_f = 2$ and $n_f = 3$ results at five and six loops using the approximate $c_{44}^{IS}$ coefficient from Eq. (4.10). Same caption as Table IV.</caption><thead><tr><th>RGOPT [$\rho(\lambda, \alpha_S)$] $\delta^k$, RG order</th><th>ln<sup>2</sup><sub>μ</sub></th><th>α<sub>S</sub></th><th>-⟨q̃q⟩<sup>1/3</sup><sub>RGI</sub> / λ<sub>2</sub></th><th>μ / λ<sub>2</sub></th></tr></thead><tbody><tr><td colspan="5">n<sub>f</sub> = 2:</td></tr><tr><td>δ<sup>4</sup>, RG four loop (full)</td><td>-0.71996 ± 0.23156i</td><td>0.37549 ± 0.05979i</td><td>0.79341 ± 0.0099350i</td><td>3.4177</td></tr><tr><td>δ<sup>4</sup>, RG five loop (full)</td><td>-0.74047 ± 0.24939i</td><td>0.36580 ± 0.053884i</td><td>0.79353 ± 0.0098374i</td><td>3.5523</td></tr><tr><td>δ<sup>4</sup>, RG four loop (α<sub>S</sub><sup>1/2</sup>)</td><td>-0.92679 ± 0.10146i</td><td>0.39942 ± 0.0074906i</td><td>0.79090 ± 0.0015884i</td><td>3.3590</td></tr><tr><td>δ<sup>4</sup>, RG five loop (α<sub>S</sub><sup>1/2</sup>)</td><td>-0.926667 ± 0.10103i</td><td>0.39949 ± 0.0074377i</td><td>0.79166 ± 0.0016397i</td><td>3.3583</td></tr><tr><td>δ<sup>5</sup>, RG five loop (full)</td><td>-1.15411</td><td>0.32270</td><td>0.75896</td><td>4.4605</td></tr><tr><td colspan="5">n<sub>f</sub> = 3:</td></tr><tr><td>δ<sup>4</sup>, RG four loop (full)</td><td>-0.77002 ± 0.19327i</td><td>0.40513 ± 0.082108i</td><td>0.75375 ± 0.0019729i</td><td>3.3255</td></tr><tr><td>δ<sup>4</sup>, RG five loop (full)</td><td>-0.79824 ± 0.18251i</td><td>0.40213 ± 0.067656i</td><td>0.75291 ± 0.00097752i</td><td>3.4381</td></tr><tr><td>δ<sup>4</sup>, RG four loop (α<sub>S</sub><sup>1/2</sup>)</td><td>-1.1570</td><td>0.42880</td><td>0.76070</td><td>3.3996</td></tr><tr><td>δ<sup>4</sup>, RG five loop (α<sub>S</sub><sup>1/2</sup>)</td><td>-1.1572</td><td>0.42885</td><td>0.75957</td><td>3.3905</td></tr><tr><td>δ<sup>5</sup>, RG five loop (full)</td><td>-1.3589</td><td>0.33097</td><td>0.71983</td><td>4.8713</td></tr></tbody></table>
<table><caption>TABLE VI. $n_f = 0$ (quenched approximation) results up to five loops; same caption as Table I.</caption><thead><tr><th>δ<sup>k</sup> order</th><th>ln(μ̃<sup>2</sup>)</th><th>α̃<sub>S</sub></th><th>-⟨⟨q̃q⟩⟩<sup>1/3</sup>/λ<sub>0</sub></th><th>μ̃/λ<sub>0</sub></th></tr></thead><tbody><tr><td>δ, RG one loop</td><td>-179/1452</td><td>11π/40 ≈ 0.86394</td><td>1.0072</td><td>1.9370</td></tr><tr><td>δ, RG two loop</td><td>-0.27548</td><td>0.49747</td><td>0.9122</td><td>2.4945</td></tr><tr><td>δ<sup>2</sup>, RG two loop</td><td>-0.54201</td><td>0.47091</td><td>0.85679</td><td>2.5795</td></tr><tr><td>δ<sup>2</sup>, RG three loop</td><td>-0.54434</td><td>0.41691</td><td>0.83522</td><td>2.8224</td></tr><tr><td>δ<sup>3</sup>, RG three loop</td><td>-0.63438</td><td>0.39314</td><td>0.83741</td><td>2.9713</td></tr><tr><td>δ<sup>3</sup>, RG four loop</td><td>-0.60347</td><td>0.38002</td><td>0.83508</td><td>3.0687</td></tr><tr><td>δ<sup>4</sup>, RG four loop</td><td>-0.54077</td><td>0.39241</td><td>0.85253</td><td>2.9764</td></tr><tr><td>δ<sup>4</sup>, RG five loop</td><td>-0.45296</td><td>0.42275</td><td>0.86069</td><td>2.8045</td></tr></tbody></table>

simulations as we will examine. Specializing our calculations to the quenched approximation from the known exact $n_f$ dependence in all relevant perturbative and RG coefficients, one simply takes $n_f = 0$ everywhere consistently. Proceeding as previously described, from the first NLO (two-loop) nontrivial order up to five loops gives the results in Table VI.

The solutions in Table VI for the quenched case, up to five-loop order included, are all located on the real AF-matching branch (which is unique at a given order), although when using five-loop RG in the very last line, the solution is located very close to the border of non-AF-matching branches. We observe that the condensate magnitude $|\langle\bar{q}q\rangle|^{1/3}$ is driven to ~2.5% higher values when going from four to five loops, as similarly observed above for $n_f = 2, 3$ in Tables III and IV. As above mentioned this is essentially traced to the instability from the first occurrence of $\pi^4$ discontinuity term at five loops.

Although in this quenched case one obtains at five loops a real solution upon using the complete RG, it is also instructive to examine the trend obtained from successive RG truncations, or alternatively when performing the calculation at six loops, giving the results in Table VII. As one can see here it is the effect of RG truncation that pushes the AF-matching solution to (slightly) nonreal values. Comparing these RG-truncated results from

Table VII, which have negligible imaginary parts, with the corresponding real solutions using the complete RG in Table VI, gives a useful estimate of the impact of such RG-consistent truncations. Concerning the six-loop results, very similar to the $n_f = 2$ and $n_f = 3$ cases they are real and very regular, and again much closer to the four-loop results in Table VI.

Finally for completeness we have also considered the $n_f = 1$ case: Although it is not very relevant physically, it can be viewed at least as a further consistency cross-check of our results. We have explored variants similar to other $n_f$ values above but simply summarize here the main results. At four-loop order one obtains the unique real solution

$$ \langle\bar{q}q\rangle_{\text{RGI}}^{1/3}(n_f = 1, \text{four loop}) = -0.8039\bar{\Lambda}_1 \quad (5.4) $$

which appears very close to the “average” of $n_f = 0$ and $n_f = 2$ four-loop results. At five loops, using four- or five-loop RG, the real solution is no longer on the AF-matching branch, similar to the $n_f = 2, 3$ cases. The unique AF-matching solution obtained from truncating $\alpha_S^{k\ge 6}$ in the RG equation gives

$$ \begin{aligned} \langle\bar{q}q\rangle_{\text{RGI}}^{1/3}(n_f = 1, \text{truncated RG five loop}) \\ = -(0.8271 \pm 0.0007i)\bar{\Lambda}_1. \end{aligned} \quad (5.5) $$

Finally at six loops a real solution is recovered, giving

$$ \langle\bar{q}q\rangle_{\text{RGI}}^{1/3}(n_f = 1, \text{RG five loop}) = -0.7984\bar{\Lambda}_1. \quad (5.6) $$

Similar to the $n_f = 0, 2, 3$ cases, once more the five-loop results produce a substantial ~2.8% increase of the condensate as compared to the four-loops results, while the six-loop results are very close to the latter.

## F. Evaluating theoretical uncertainties

Comparing the different above results from $n_f = 0$ to $n_f = 3$, it is tempting to consider the manifestly more stable results obtained at four loops and six loops as a likely better approximation than the more sensibly shifted

<table><caption>TABLE VII. $n_f = 0$ results at five and six loops; same caption as Table III.</caption><thead><tr><th>RGOPT [ρ(λ, α<sub>S</sub>)] δ<sup>k</sup>, RG order</th><th>ln(μ̃<sup>2</sup>)</th><th>α̃<sub>S</sub></th><th>-⟨⟨q̃q⟩⟩<sup>1/3</sup>/λ<sub>0</sub></th><th>μ̃/λ<sub>0</sub></th></tr></thead><tbody><tr><td>δ<sup>4</sup>, RG four loop (full)</td><td>-0.54077</td><td>0.39241</td><td>0.85253</td><td>2.9764</td></tr><tr><td>δ<sup>4</sup>, RG five loop (full)</td><td>-0.45296</td><td>0.42275</td><td>0.86069</td><td>2.8045</td></tr><tr><td>δ<sup>4</sup>, RG four loop (α<sub>S</sub><sup>k≥6</sup>)</td><td>-0.70122 ± 0.11373i</td><td>0.36891 ± 0.0031072i</td><td>0.84724 ± 0.00030466i</td><td>3.1610</td></tr><tr><td>δ<sup>4</sup>, RG five loop (α<sub>S</sub><sup>k≥6</sup>)</td><td>-0.70089 ± 0.11265i</td><td>0.36906 ± 0.0030291i</td><td>0.85101 ± 0.00019590i</td><td>3.170</td></tr><tr><td>δ<sup>5</sup>, RG five loop (full)</td><td>-0.90052</td><td>0.30873</td><td>0.82242</td><td>3.912</td></tr><tr><td>δ<sup>5</sup>, RG five loop (α<sub>S</sub><sup>k≥6</sup>)</td><td>-0.97183</td><td>0.31740</td><td>0.82383</td><td>3.772</td></tr></tbody></table>
five-loop results. Note indeed that if discarding the latter, the combined four-loop and six-loop results would provide a seemingly very accurate determination. But since the six-loop results are only partial, we more conservatively combine all those results within our estimate of uncertainties. More precisely we take the average between the four-, five- and six-loop results as our central values and their differences as our theoretical uncertainties (taking at five loops the real parts of the results having the smallest imaginary parts, that are consequently more reliable). Then we estimate the uncertainties linearly from the complete range spanned by maximal and minimal values.

For $n_f = 0$, for which real solutions occur at all RGOPT successive orders considered, we obtain

$$ \langle \bar{q}q \rangle_{\text{RGI}}^{1/3}(n_f = 0) \approx -(0.840^{+0.020}_{-0.016}) \bar{\Lambda}_0, \quad (5.7) $$

where we give only a three digit accuracy given the uncertainties.

For $n_f = 2$, proceeding similarly we obtain

$$ \langle \bar{q}q \rangle_{\text{RGI}}^{1/3}(n_f = 2) \approx -(0.781^{+0.019}_{-0.009}) \bar{\Lambda}_2. \quad (5.8) $$

And finally for $n_f = 3$:

$$ \langle \bar{q}q \rangle_{\text{RGI}}^{1/3}(n_f = 3) \approx -(0.751^{+0.019}_{-0.010}) \bar{\Lambda}_3. \quad (5.9) $$

Equations (5.7)–(5.9) constitute our primary results, as these do not depend on any extra theoretical or experimental input besides the basic perturbative content used in the calculation and the RGOPT method. Now, to make contact with other independent determinations of the quark condensate, often conventionally given at the standard scale $\mu \approx 2$ GeV for reference, one needs to perform a (perturbative) renormalization scale evolution. One should keep in mind that, in contrast with the above results, such RG evolution unavoidably also entails $\alpha_s$ (and other related) uncertainties.

## VI. $\langle \bar{q}q \rangle (\mu = 2 \text{ GeV})$ AND COMPARISON WITH OTHER DETERMINATIONS

To evolve perturbatively the condensate from our results above, the simplest procedure is to take the values obtained for the scale-invariant condensate (5.3), within uncertainties, Eqs. (5.7)–(5.9) and extract from these the condensate at another chosen (perturbative) scale $\mu'$, using again (5.3) at five-loop order, now taking $g \equiv 4\pi\alpha_s(\mu')$, after evolving $\alpha_s(\mu)$ at five-loop order of Eq. (5.2) toward the conventional scale $\mu' = 2$ GeV. The overall reliability of this (perturbative) evolution is to be assessed on the grounds that the primary RGOPT results above at four, five and six loops are obtained at reasonably perturbative optimized scale values ($3.3\Lambda_{\overline{\text{MS}}} \lesssim \tilde{\mu} \lesssim 4.8\Lambda_{\overline{\text{MS}}}$; compare Tables I–IV).

It is more appropriate to separate the discussion below for different $n_f$ values, since those do not have all the same reliability status (also when comparing our results with other independent determinations of the condensate) as we discuss next. We consider successively $n_f = 3$, $n_f = 2$, and $n_f = 0$ (quenched approximation).

### A. $n_f = 3$

For $n_f = 3$, one can use very reliable $\alpha_s$ determinations in the perturbative range. We also account properly for the charm quark mass threshold effects [53] on $\alpha_s(\mu \sim m_c)$. From the most recent world average $\alpha_s(m_Z)$ value [51]

$$ \alpha_s(m_Z) = 0.1179 \pm 0.0010, \quad (6.1) $$

we obtain in a first stage, accounting for threshold effects at $\mu \sim m_b$ and $\mu \sim m_c$,$^{11}$

$$ \bar{\Lambda}(n_f = 3) = (331 \pm 16) \text{ MeV} \quad (6.2) $$

and

$$ \alpha_s(2 \text{ GeV}) = 0.3007 \pm 0.008. \quad (6.3) $$

Then using Eq. (5.3) applied to Eq. (5.9) leads to

$$ \langle \bar{q}q \rangle_{n_f=3}^{1/3}(2 \text{ GeV}) = -(0.826^{+0.021}_{-0.011}) \bar{\Lambda}_3. \quad (6.4) $$

Thus combining Eq. (6.2) with (6.4) leads to

$$ \langle \bar{q}q \rangle_{n_f=3}^{1/3}(2 \text{ GeV}, \bar{\Lambda}_3^{\text{wa}}) \approx -(273^{+7}_{-4} \pm 13) \text{ MeV} \quad (6.5) $$

where the first error is our rather conservative theoretical RGOPT uncertainty from Eq. (6.4) and the second one is from $\bar{\Lambda}(3)$ uncertainty. (Since these two uncertainties have very different origin we do not combine them.) It is worth remarking at this point that Eq. (6.4) is only slightly shifted with respect to our previous (average of three- and four-loop RGOPT) result [26], while the central value and uncertainties in Eq. (6.5) [compare Eq. (6.6) of [26]] are principally affected by the slight decrease of the most recent $\alpha_s$ world average with substantial increase of uncertainties (see [51] for detailed explanations on these features).

<footnote>$^{11}$For the rather low values of the scales involved, it appears more appropriate to use the exponentiated forms Eqs. (5.2) and (5.3), somewhat more stable than their purely perturbative expansions. We have also cross-checked our five-loop RG evolution with the results using the well-known public code RUNDEC [54] recently upgraded to five-loop order.</footnote>
To compare with other independent determinations, first the most precise $n_f = 3$ lattice determination we are aware of, in the chiral limit, is $\langle \bar{q}q \rangle_{n_f=3}^{1/3}(2 \text{ GeV}) = -(245 \pm 5 \pm 8) \text{ MeV}$ [55]. Our results are thus marginally compatible with the latter, within uncertainties of both results. Note, however, that various recent lattice results vary in a wider range for $n_f = 3$, as compiled from [2]: from $214 \pm 6 \pm 24$ [56] to $290 \pm 15$ [57]. This is largely due to the still difficult required extrapolation of lattice results to the chiral limit, which for the SU(3) case is affected by large uncertainties. A recent very precise $n_f = 3$ lattice calculation [58] using time moments of a heavy-strange pseudoscalar correlator, has obtained $\langle \bar{s}s \rangle^{1/3}(2 \text{ GeV}) \sim -(296 \pm 11) \text{ MeV}$. Since it is not in the chiral limit, it should not be directly compared with our result, given the large strange quark mass involved. Indeed as our $n_f = 3$ results are based on a relatively accurate RGOPT determination (5.9), and (6.5) obtained from a reliable $\bar{\Lambda}(3)$ world average, they appear useful independent determinations since being in the strict chiral limit, thus relevant to possibly assess the actual impact from explicit chiral symmetry breaking by the strange quark mass, by comparison with other determinations that include the latter, like [58].

## B. $n_f = 2$

For $n_f = 2$, one cannot directly link our results to the true phenomenological perturbative range values of $\alpha_S$ as above. Nevertheless, given that the (optimized) coupling values obtained in Tables I and III are reasonably perturbative, we can consider a perturbative (five-loop) RG evolution (consistently performed in a simplified QCD picture where the strange and heavier quarks are all infinitely massive, i.e., “integrated out”). To give a final (numerical) determination of the condensate, we need a value for $\bar{\Lambda}(n_f = 2)$. To our knowledge there are not many nonperturbative results for $\bar{\Lambda}(2)$ (as compared with the numerous studies for $n_f = 3$), and those results mostly originate from lattice calculations. We therefore rely on a lattice determination [59] (that best fulfill the reliability criteria of the review [2]) obtained from the Schrödinger functional method:

$$ \bar{\Lambda}(n_f = 2) = (310 \pm 20) \text{ MeV}. \quad (6.6) $$

One should keep in mind, however, that somewhat larger uncertainties are obtained if taking more conservatively all presently available lattice results [59,60], as compiled in [2].<fnref n="12" /> After RG evolution up to 2 GeV we obtain accordingly from Eq. (5.8):

$$ \langle \bar{q}q \rangle_{n_f=2}^{1/3}(2 \text{ GeV}) = -(0.863^{+0.021}_{-0.010}) \bar{\Lambda}_2. \quad (6.7) $$

[Notice that this strictly $n_f = 2$ result, thus with the strange and heavier quarks integrated out, correspondingly has $\alpha_S(\mu)$ values not consistent with the phenomenological values of Eq. (6.3): Instead we find $\alpha_S^{n_f=2}(2 \text{ GeV}) \approx 0.262 \pm .008$. Combining (6.7) with Eq. (6.6) leads to<fnref n="13" />

$$ \langle \bar{q}q \rangle_{n_f=2}^{1/3}(2 \text{ GeV}, \bar{\Lambda}_2) \approx -(267^{+7}_{-4} \pm 18) \text{ MeV}, \quad (6.8) $$

where again the first error range is our RGOPT uncertainty from Eq. (6.7) while the second one is from the $\bar{\Lambda}(2)$ uncertainty. Since lattice uncertainties are mostly statistical and systematic, while ours are theoretical, it is not obvious to combine these in a sensible manner and we keep more conservatively separate uncertainties.

To compare our result (6.8) with other recent determinations, first the presumably most precise $n_f = 2$ lattice determination to date is also from the spectral density [15] $\langle \bar{q}q \rangle_{n_f=2}^{1/3}(\mu = 2 \text{ GeV}) = -(261 \pm 6 \pm 8)$, where the first error is statistical and the second is systematic. Our results are thus very compatible within uncertainties. Note, however, that the above quoted lattice value [15] was obtained by fixing the scale with the kaon decay constant $F_K$ determined in the quenched approximation. Overall, recent $n_f = 2$ lattice determinations of the condensate in the chiral limit from several independent methods are much more precise than those for $n_f = 3$. We quote the estimate recently performed in [2], by combining results from [15,62]: $|\langle \bar{q}q \rangle|_{n_f=2}^{1/3} = -(266 \pm 10) \text{ MeV}$, where the uncertainties include both systematic and statistical ones.

One may also compare with recent results from spectral sum rules [4]: $\langle \bar{u}u \rangle^{1/3}(\equiv \langle \bar{d}d \rangle^{1/3})(2 \text{ GeV}) \approx -(276 \pm 7) \text{ MeV}$, but keeping in mind that the latter sum rules actually determine precisely the current quark masses, so that the $\langle \bar{u}u \rangle$ value is indirectly extracted from using the GMOR relation (1.1). Accordingly the comparison is not strictly for the chiral limit. In this context, even though the overall reliability of $\bar{\Lambda}(2)$ is not yet at the level of $\bar{\Lambda}(3)$, the results (5.8) and (6.8) constitute reasonably accurate independent determinations in the chiral limit. Indeed, given that the present $n_f = 2$ lattice results for the condensate are quite accurate, it is tempting alternatively to combine the latter with our firmer result Eq. (5.8), in order to rather determine a new independent estimate of $\bar{\Lambda}(2)$:

<footnote n="12">¹²Our own determination of $\bar{\Lambda}(n_f = 2)$ from the pion decay constant $F_\pi$ using three-loop RGOPT [25] is compatible with the range obtained from lattice simulations, but also has larger uncertainties, as compared to our RGOPT $\bar{\Lambda}(n_f = 3)$ results.</footnote>

<footnote n="13">¹³As compared with our 2015 four-loop result [26], note that the central value in Eq. (6.8) is principally affected by the somewhat lower central $\bar{\Lambda}(n_f = 2)$ value from (6.6) as compared with previously used value from [61].</footnote>
Taking the above quoted estimate of the condensate given by [2], this gives

$$ \bar{\Lambda}(n_f = 2) = (308_{-6}^{+4} \pm 12) \text{ MeV.} \quad (6.9) $$

### C. $n_f = 0$ (quenched approximation)

Finally for completeness we also give results for the quenched approximation ($n_f = 0$). In this case we evolve $\alpha_s(\mu)$ and use Eq. (5.3) at five loops but in the appropriate $n_f = 0$ approximation. As previously we need $\bar{\Lambda}(n_f = 0)$, available from various approaches with lattice simulations. We rely on the average performed in [2], combining different precise lattice results [63]:

$$ \bar{\Lambda}(n_f = 0) = (257 \pm 7) \text{ MeV.} \quad (6.10) $$

As stressed in [2], it is worth noting that this value is obtained by using the same value as for $n_f = 2$ and $n_f = 2 + 1$ of the basic lattice scale defined from the quark static potential, $r_0 = 0.472$ fm, which for $n_f = 0$ amounts merely to a defining convention for $\bar{\Lambda}(n_f = 0)$.

Next the RG evolution from Eq. (6.12) leads to

$$ \langle \bar{q}q \rangle_{n_f=0}^{1/3} (2 \text{ GeV}) = -(0.932_{-0.018}^{+0.022}) \bar{\Lambda}_0. \quad (6.11) $$

Combining (6.10) with Eq. (6.11) we obtain

$$ \langle \bar{q}q \rangle_{n_f=0}^{1/3} (2 \text{ GeV}, \bar{\Lambda}_0) = -(240_{-5}^{+6} \pm 6) \text{ MeV.} \quad (6.12) $$

It appears to us not easy to compare (6.12) with other determinations, since most phenomenological determinations of the condensate are obviously obtained for $n_f \ge 2$. Concerning lattice simulations, most of the modern calculations no longer use the quenched approximation, performing simulations with fully dynamical sea quarks, while too old results in the quenched approximation are presumably affected by rather large uncertainties. To our knowledge, there is one precise, often quoted latest quenched simulation result [64]:

$$ \langle \bar{q}q \rangle_{n_f=0}^{1/3} (2 \text{ GeV, lattice}) = -(250 \pm 3) \text{ MeV.} \quad (6.13) $$

So our result (6.12) appears consistent with the latter within uncertainties. We stress, however, that our study of the quenched case $n_f = 0$ is merely motivated as a consistency cross-check of our method, since the quenched approximation is anyway not very realistic.

### D. Further discussion on $\langle \bar{q}q \rangle(n_f)$ dependence

Comparing all our results for $n_f = 0, 2, 3$ at the same perturbative orders, it appears that the ratio of the quark condensate to $\bar{\Lambda}^3$ has a sizable but moderate dependence on the number of flavors $n_f$: There is a clear trend that $|\langle \bar{q}q \rangle^{1/3}|(n_f)/\bar{\Lambda}(n_f)$ decreases regularly, roughly linearly by about 4% for $n_f \to n_f + 1$ (that is clear at least from the studied $n_f \le 3$ cases). Naively (perturbatively) the moderate dependence on $n_f$ is expected (as long as $n_f$ is not large), since it only appears explicitly at three-loop order. Nevertheless it does not imply a similar decrease of the absolute condensate values, as those depend on $\bar{\Lambda}(n_f)$, that appears rather to increase with $n_f$ for $n_f \le 3$ [at least if considering the low lattice $\bar{\Lambda}(0)$ value (6.10), but it is not so clear from comparing $\bar{\Lambda}(2)$ and $\bar{\Lambda}(3)$ given all the present uncertainties in their values]. Concerning the $n_f = 3$ to $n_f = 2$ condensate ratio, various lattice results have still rather large uncertainties at present [2] but some recent results are more compatible with a ratio unity [58,65]. The spectral sum rules prediction for the ratio is also not very precise [66,67]: $\langle \bar{s}s \rangle / \langle \bar{u}u \rangle = 0.74_{-0.12}^{+0.34}$, (see also the recent review [68]). Since our results are by construction valid in the strict chiral limit, taken at face value they indicate that the possibly larger difference obtained by some other determinations [2,22] is more likely due to the explicit breaking from the large strange quark mass, rather than an intrinsically strong $n_f$ dependence of the condensate in the exact chiral limit.

## VII. SUMMARY AND CONCLUSION

We have reconsidered our variational RGOPT approach applied to the spectral density of the Dirac operator, the latter being obtained in a first stage from the perturbative logarithmic discontinuities of the quark condensate in the $\overline{MS}$ scheme. This construction allows successive sequences of nontrivial variationally optimized results in the strict chiral limit, from two- to five-loop levels using exactly known perturbative content, and partially up to six loops, the latter more approximately relying on the six-loop content exactly predictable from five-loop renormalization group properties. The results Eqs. (5.7)–(5.9) are those that we consider the firmer, while latter results in Eqs. (6.8) and (6.5) are further affected by present uncertainties in perturbative evolution and $\bar{\Lambda}$ values. Equations (5.7)–(5.9) show a very good stability and empirical convergence, although the strictly five-loop results exhibit some instabilities with respect to both four- and six-loop results. Those instabilities are traced to specific features of the spectral density, namely the occurrence at growing orders of new large $\pi^{2k}$ discontinuity contributions that tend to destabilize the original perturbative coefficients when first appearing at a given order. For all the considered cases from $n_f = 0$ (quenched approximation) to $n_f = 3$, it is striking that the six-loop
results are very close to the four-loop ones, both exhibiting
very stable properties. It appears convincing to us that the
systematically ~2% higher values of |⟨q̄q⟩|<sup>1/3</sup><sub>RGI</sub>(n<sub>f</sub>) obtained
at five-loop RGOPT order are largely an artifact of the
instability from the π<sup>4</sup> discontinuity terms appearing first at
five-loop order. Nevertheless we have incorporated more
conservatively the differences between the four-, five- and
six-loop results as intrinsic theoretical uncertainties, which
are of order ±2%. Notice that, if less conservatively
discarding the presumably less reliable strictly five-loop
results from our averages, the lowest values in Eqs. (5.7)–
(5.9) are favored, with much smaller uncertainties. In any
case the final condensate values and uncertainties in
Eqs. (6.5) and (6.8) are more affected by the present
uncertainties on the basic QCD scale Λ, both for n<sub>f</sub> = 2
and n<sub>f</sub> = 3. (To possibly get rid of Λ uncertainties, par-
ticularly for n<sub>f</sub> = 2, one could in principle apply RGOPT
directly to a more physical RG-invariant quantity, like
⟨q̄q⟩<sub>RGI</sub>/F<sup>3</sup><sub>π</sub>, combining the present analysis with the one
in [25] for F<sub>π</sub>, but this involves somewhat nontrivial issues
and is left for future investigation.)

In conclusion the chiral condensate values obtained in
our analysis are very compatible, within uncertainties, with
the most precise recent lattice determinations for all
considered n<sub>f</sub> values. Our results for n<sub>f</sub> = 3 are perhaps
of particular interest, given that other independent deter-
minations are either not in the chiral limit or, concerning
lattice results, are affected by still rather important uncer-
tainties in the chiral extrapolation [2]. Finally our results
have indicated a moderate flavor dependence of the
⟨q̄q⟩<sup>1/3</sup><sub>n<sub>f</sub></sub>/Λ<sub>n<sub>f</sub></sub> values in the chiral limit for n<sub>f</sub> ≤ 3.

ACKNOWLEDGMENTS

We are very grateful to Konstantin Chetyrkin for
valuable exchanges, and to him and Pavel Baikov for
providing us their results of the five-loop vacuum anoma-
lous dimension [32] before it was published. We are also
very grateful to Andreas Maier and Peter Marquard for
valuable discussions at the RADCOR 2019 conference, and
for providing us four-loop results [45] leading to Eq. (4.11)
in our normalization.

APPENDIX: RG AND OTHER PERTURBATIVE
QUANTITIES

In this Appendix we give for completeness all the
relevant quantities related to perturbative RG properties
used in our calculations. The RG coefficients up to five
loops for a general gauge theory were obtained in complete
analytical form respectively in Ref. [30] for the beta
function and [31] for the anomalous mass dimension.
We do not repeat those expressions explicitly here, refer-
ring to these articles. Note simply that we mainly use the
normalization g ≡ 4πα_s(μ) ≡ 4π^2a_s(μ), such that

$$
\beta(g) \equiv \frac{dg}{d \ln \mu} = -2 \sum_k b_k g^{k+2}, \qquad \gamma_m(g) = \sum_k \gamma_k g^{k+1}, \tag{A1}
$$

where our $b_k$, $\gamma_k$ expressions are related as $b_k = \beta_k(4\pi^2)^{-1-k}$, $\gamma_k = 2(4\pi^2)^{-1-k}\gamma_k$ with respect e.g., to the first Refs. [30,31] respectively.

1. Vacuum-energy anomalous dimension

Next, in our normalization conventions the anomalous
dimension of the vacuum energy, Γ⁰(a_s) entering Eq. (4.2),
is given to five-loop order as

$$
\Gamma^0(a_s) = -2 \times \left( \frac{3}{16\pi^2} \right) \left[ 1 + \sum_{k=1}^{4} a_s^k \Gamma_k^0 + \mathcal{O}(a_s^5) \right] \quad (\text{A2})
$$

with the coefficients up to three loops determined long
ago [42],

$$
\Gamma_1^0 = \frac{4}{3}, \qquad \Gamma_2^0 = \frac{457}{72} - \frac{5}{12}n_l - \frac{29}{12}n_h - \frac{2}{3}z_3, \qquad (\text{A3})
$$

and the four-loop Γ<sup>0</sup><sub>3</sub> and five-loop Γ<sup>0</sup><sub>4</sub> coefficients obtained
in full analytical form in Ref. [32] are given explicitly
respectively in Eqs. (3.4) and (3.5) of [32], for n<sub>f</sub> ≡ n<sub>l</sub> +
n<sub>h</sub> quark flavors and n<sub>h</sub> = 1 heavy quark.<sup>14</sup> For complete-
ness we give here their relevant expressions adapted to our
case (in numerical form for short), where in practice we
consider n<sub>l</sub> = 0 and n<sub>h</sub> (≡n<sub>f</sub>) massive degenerate quarks:

$$
\begin{equation}
\begin{split}
\Gamma_3^0(n_l, n_h) \approx {}& 33.6625 + 0.18139n_h^2 - 32.5586n_h \\
& + 0.214632n_hn_l - 4.96507n_l + 0.0332417n_l^2,
\end{split}
\tag{A4}
\end{equation}
$$

and

$$
\begin{equation}
\begin{split}
\Gamma_4^0(n_l, n_h) \approx {}& 242.021 + 0.0185562n_h^3 + 17.6037n_h^2 \\
& + 0.03715n_h^2n_l - 43.3192n_l + 1.0631n_l^2 \\
& + 0.0000376492n_l^3 - 299.998n_h \\
& + 18.6668n_hn_l + 0.0186315n_hn_l^2.
\end{split}
\tag{A5}
\end{equation}
$$

2. Perturbative condensate and spectral density

Next, the coefficients of the perturbative quark condensate, Eq. (4.1), were given at three loops in Eq. (4.4). Using Eqs. (4.2) with Eq. (3.4) we determine the relevant coefficients at four-loop and higher orders:

<footnote n="14"><sup>14</sup>Note a trivial factor 2 normalization difference in Eq. (A2) with respect to [32] due to our use of d ln µ in Eq. (4.2). Also Γ<sup>0</sup><sub>k</sub> in Eqs. (A2)–(A5) differ by an overall factor 3/(16π<sup>2</sup>) from the original (γ<sup>di.nd</sup><sub>0</sub>)<sub>k</sub> in the notations of [32].</footnote>
$$c_{40} = \frac{1}{54}(81 - 2n_f)(57 - 2n_f),$$

$$c_{41} = -\frac{1}{81}(18171 + 20n_f^2 - 1353n_f),$$

$$c_{42} = \frac{1}{8 \times 81} [226647 - 30690n_h - 18378n_l + 632n_h^2 + 832n_hn_l + 200n_l^2 - 72(57 + 28n_f)z_3],$$

$$c_{43} = \frac{1}{32 \times 729} [4n_h^2(4752z_3 - 5935) - 4n_l^2(3024z_3 + 751) + 8n_l n_h(864z_3 - 3343) \\ - 12n_f(2304a_4 + 96l_{z2} - 5040z_4) + 12n_l(22428z_3 + 53093) - 12n_h(30060z_3 + 9720z_5 - 124373) \\ + 9(87552a_4 + 3648l_{z2} + 10368z_3 - 68400z_4 - 81840z_5 - 854633)], \qquad (A6)$$

where $n_f = n_l + n_h$, $l_{z2} = \ln^2(2)(\ln^2(2) - 6z_2)$, $z_i = \zeta(i)$.

In numerical approximation this gives

$$c_{40} = 85.5 - 5.11111n_f + 0.0740741n_f^2,$$

$$c_{41} = -224.333 + 16.7037n_f - 0.246914n_f^2,$$

$$c_{42} = 342.151 - 51.1008n_h - 32.1008n_l + 0.975309n_h^2 + 0.308642n_l^2 + 1.28395n_hn_l,$$

$$c_{43} = -375.082 + 42.6214n_h + 43.5949n_l - 0.0382074n_h^2 - 0.790268n_hn_l - 0.752061n_l^2, \qquad (A7)$$

while the complete expression for the last nonlogarithmic four-loop coefficient $c_{44}$ is given explicitly in Eq. (4.11). Similarly we obtain for the five-loop logarithmic coefficients:

$$c_{50} = -\frac{1}{4 \times 135}(81 - 2n_f)(57 - 2n_f)(49 - 2n_f), \qquad (A8)$$

$$c_{51} = -\frac{1}{8 \times 81}(-952101 + 40n_f^3 - 3898n_f^2 + 112749n_f), \qquad (A9)$$

$$c_{52} = \frac{1}{16 \times 243} (2150625n_h - 95963n_h^2 + 1264n_h^3 - 146134n_hn_l \\ + 2928n_h^2n_l - 50171n_l^2 + 2064n_hn_l^2 + 400n_l^3 \\ - 144n_f(-869 + 28n_f)z_3 + 3(515779n_l + 72(-56554 + 931z_3))), \qquad (A10)$$

$$729 \times 4^3 c_{53} = 32n_h^3(648z_3 - 773) - 32n_l^3(324z_3 + 125) \\ + 9(-768a_4 + 32l_{z2})(57 - 2n_f)(49 - 2n_f) \\ + 116448z_3 + 1675800z_4 + 1292280z_5 + 23864201) \\ + n_h(-4320(27n_f - 836)z_5 + 6310872z_3 - 1455840z_4 - 50009817) \\ + 96n_l n_h^2(324z_3 - 557) + n_l n_h(-131328z_3 + 43200z_4 + 3132538) \\ - 15n_l(608040z_3 + 97056z_4 - 50256z_5 + 1871063) \\ + n_h^2(-761616z_3 + 21600z_4 + 2316653) \\ + n_l^2(-32736n_h + 630288z_3 + 21600z_4 + 815885), \qquad (A11)$$
$$ 
\begin{aligned}
243 \times 2 \times 4^3 c_{54} = & -243 \times 4^3 (49 - 2n_f)c_{44} \\
& + n_h^3(48z_3 - 1728z_4 + 875) + n_l^3(48z_3 - 1728z_4 + 1451) \\
& + n_h^2 n_l(144z_3 - 5184z_4 + 3201) \\
& + n_h n_l^2(144z_3 - 5184z_4 + 3777) \\
& + n_h n_l \left( -144z_3(81z_3 + 2735) + 128736z_4 + 268992z_5 - 48600z_6 - \frac{9903989}{18} \right) \\
& + n_h^2 \left( -72z_3(162z_3 + 3191) + 96444z_4 + 173376z_5 - 48600z_6 - \frac{16643189}{36} \right) \\
& + n_l^2 \left( -164088z_3 + 32292z_4 + 95616z_5 - \frac{3164789}{36} \right) \\
& + n_h \left( -162816a_4 - 6784l_{z2} + 12z_3(39228z_3 + 407969) - 932766z_4 - 8087688z_5 + 262800z_6 + 3213000z_7 + \frac{54045443}{6} \right) \\
& + n_l \left( -162816a_4 - 6784l_{z2} + 12(134429 - 12990z_3)z_3 + 638958z_4 - 461448z_5 - 927900z_6 - 95256z_7 + \frac{21934883}{6} \right) \\
& + 1170432a_4 + 48768l_{z2} + 12z_3(8031z_3 + 1007870) - 5759550z_4 - 27852768z_5 + 10192050z_6 + 8641836z_7 - \frac{145813179}{8},
\end{aligned}
\quad (A12) $$

where we conveniently expressed $c_{54}$ in terms of the four-loop nonlogarithmic $c_{44}$ coefficient. In numerical approximation we obtain

$$ c_{50} = -418.95 + 42.1444n_f - 1.38519n_f^2 + 0.0148148n_f^3, $$

$$ c_{51} = 1469.29 - 173.995n_f + 6.01543n_f^2 - 0.0617284n_f^3, $$

$$ c_{52} = -3079.72 + 436.666n_f - 14.1506n_f^2 + 0.102881n_f^3 \\ + n_h(155.167 - 11.7778n_f + 0.222222n_f^2), $$

$$ c_{53} = 5102.45 - 852.446n_h - 843.205n_l + 61.7769n_hn_l \\ + 27.7449n_h^2 + 34.032n_l^2 - 0.344719n_h^2n_l \\ - 0.701646n_hn_l^2 + 0.00406919n_h^3 - 0.352858n_l^3, $$

$$ c_{54} = (n_f - 24.50)c_{44} - 617.146 + 309.613n_h + 144.324n_l \\ - 16.7381n_h^2 - 4.8565n_l^2 - 21.5946n_hn_l \\ - 0.0719093n_h^2n_l - 0.0533908n_hn_l^2 - 0.0301426n_h^3 - 0.0116241n_l^3. \quad (A13) $$

It is straightforward to apply the RG Eq. (3.4) to obtain similarly the six-loop logarithmic coefficients $c_{6k}$ of $L_m^{6-k}$. To avoid unnecessarily lengthy expressions it is convenient to equivalently express the $c_{6k}$ as functions of the above lower-order perturbative and RG coefficients: With the same normalization as in Eq. (4.1), with an overall $3/(2\pi^2)m^3a_s^5$ factor at six loops, they read
$$
\begin{align}
c_{60} &= -\frac{8\pi^2}{3}(2b_0 + \gamma_0), \nonumber \\
c_{61} &= -\frac{4\pi^2}{5}(4(2b_0 + \gamma_0)c_{51} + 5c_{50}\gamma_0 + 8\pi^2 c_{40}(3b_1 + 2\gamma_1)), \nonumber \\
c_{62} &= -4\pi^2[(2b_0 + \gamma_0)c_{52} + \gamma_0 c_{51} - 2\pi^2((3b_1 + 2\gamma_1)c_{41} + 2c_{40}\gamma_1) \nonumber \\
&\quad - 16\pi^4 c_{30}(b_2 + \gamma_2)], \nonumber \\
c_{63} &= -\frac{4}{3}\pi^2 \Biggl[ 3\gamma_0 c_{52} + 4(2b_0 + \gamma_0)c_{53} + 4\pi^2(2(3b_1 + 2\gamma_1)c_{42} + 3\gamma_1 c_{41} \nonumber \\
&\quad + 4\pi^2(3c_{30}\gamma_2 + 4c_{31}(b_2 + \gamma_2) + 8\pi^2 c_{20}(b_3 + 2\gamma_3))) \Biggr], \nonumber \\
c_{64} &= 4\pi^2[-2(2b_0 + \gamma_0)c_{54} - \gamma_0 c_{53} - 4\pi^2((3b_1 + 2\gamma_1)c_{43} + \gamma_1 c_{42}) \nonumber \\
&\quad - 16\pi^4(2c_{32}(b_2 + \gamma_2) + c_{31}\gamma_2) \nonumber \\
&\quad - 64\pi^6((b_3 + 2\gamma_3)c_{21} + \gamma_3 c_{20}) + 512\pi^8\gamma_4]. \tag{A14}
\end{align}
$$

Next from Eqs. (4.15) and (4.16) it is straightforward to derive the corresponding six-loop coefficients of the spectral density in the normalization of Eq. (4.17), that we give here for completeness:

$$
\begin{align}
\rho_{61} &= 3c_{60}, & \rho_{62} &= \frac{5}{2}c_{61}, & \rho_{63} &= 2c_{62} - \frac{5\pi^2}{2}c_{60}, & \rho_{64} &= \frac{3}{2}c_{63} - \frac{5\pi^2}{4}c_{61}, \nonumber \\
\rho_{65} &= c_{64} - \frac{\pi^2}{2}c_{62} + \frac{3\pi^4}{16}c_{60}, & \rho_{66} &= \frac{1}{2}c_{65} - \frac{\pi^2}{8}c_{63} + \frac{\pi^4}{32}c_{61}. \tag{A15}
\end{align}
$$

Note that from standard RG properties, only $c_{65}$ [and therefore only the nonlogarithmic coefficient $\rho_{66}$ of $\rho(\lambda)$ in Eq. (A15)] depends on the presently unknown six-loop vacuum energy $\Gamma_5^0$ and nonlogarithmic five-loop coefficient $c_{55}$. Accordingly, as explained in the main text, we simply ignore $c_{65}$ and $\rho_{66}$ in our six-loop analysis.

### 3. RG-invariant perturbative subtraction

Next, we also derive for completeness the coefficients entering the subtraction function $S(m, g)$ defined in Eqs. (4.12) and (4.13), such that $m\langle\bar{q}q\rangle - S(m, g)$ defines a (finite) condensate obeying the homogenous RG equation Eq. (3.4) up to five loops, that can be useful for different purposes. We obtain after some algebra<fnref n="15" />

$$
\begin{align}
s_0 &= \frac{1}{8\pi^2(b_0 - 2\gamma_0)} = -\frac{6}{(15 + 2n_f)}, \nonumber \\
s_1 &= (2\gamma_0)^{-1} \left( 4\pi^2(b_1 - 2\gamma_1) - \frac{\Gamma_1^0}{8\pi^2} \right) = \frac{633 + n_f}{48(15 + 2n_f)}, \nonumber \\
s_2 &= (2\gamma_0 + b_0)^{-1} \left( (4\pi^2)^2(b_2 - 2\gamma_2)s_0 - 2(4\pi^2)\gamma_1 s_1 - \frac{\Gamma_2^0}{8\pi^2} \right) \nonumber \\
    &= \frac{n_h(-24519 + 33408z_3) + n_l(1401 + 33408z_3) - 746n_l n_h - 2101n_h^2 + 1355n_l^2 - 27(4151 + 320z_3)}{144(15 + 2nf)(-81 + 2n_f)}, \nonumber \\
s_3 &= (2\gamma_0 + 2b_0)^{-1} \left( (4\pi^2)^3(b_3 - 2\gamma_3)s_0 - 2(4\pi^2)^2\gamma_2 s_1 - 4\pi^2(2\gamma_1 + b_1)s_2 - \frac{\Gamma_3^0}{8\pi^2} \right), \nonumber \\
s_4 &= (2\gamma_0 + 3b_0)^{-1} \left( (4\pi^2)^4(b_4 - 2\gamma_4)s_0 - 2(4\pi^2)^3\gamma_3 s_1 - (4\pi^2)^2(2\gamma_2 + b_2)s_2 - 2(4\pi^2)(\gamma_1 + b_1)s_3 - \frac{\Gamma_4^0}{8\pi^2} \right) \tag{A16}
\end{align}
$$

where we give both generic compact expressions and their particular QCD values, the latter for the four-loop $s_3$ and five-loop $s_4$ coefficients given numerically to $10^{-6}$ accuracy as

<footnote n="15">¹⁵The normalization of $s_i$ coefficients in Eq. (4.13) is different from the one in Ref. [26] for convenience, but Eq. (A16) is consistent with our previous expressions up to three-loop order in [26].</footnote>
$$s_3 = \frac{32}{(-81 + 2n_f)(-57 + 2n_f)(15 + 2n_f)} \times (2485.78 + 1045.17n_h - 1314.62n_l + 436.946n_hn_l + 351.286n_h^2 + 85.66n_l^2 - 8.6872n_h^3 - 2.56879n_l^3 - 19.9432n_h^2n_l - 13.8248n_hn_l^2 + 0.00262828n_h^2n_l^2 + 0.0573077n_h^3n_l - 0.0538034n_hn_l^3 + 0.0282158n_h^4 - 0.0273397n_l^4)}, \quad (A17)$$

$$s_4 = \frac{16}{(-49 + 2n_f)(-81 + 2n_f)(-57 + 2n_f)(15 + 2n_f)} \times (-1947880 + 0.0102918n_h^6 + 13.0258n_h^5 + 0.0432323n_h^5n_l + 682641.n_l - 98179.9n_l^2 + 3437.72n_l^3 + 42.1262n_l^4 - 1.59314n_l^5 - 0.00822672n_l^6 - 943.625n_h^4 + 50.5099n_h^4n_l + 0.0617844n_h^4n_l^2 + 19279.n_h^3 - 2788.75n_h^3n_l + 71.782n_h^3n_l^2 + 0.0206508n_h^3n_l^3 - 105289.n_h^2 + 41995.6n_h^2n_l - 2704.5n_h^2n_l^2 + 42.5442n_h^2n_l^3 - 0.0308082n_h^2n_l^4 - 723426.n_h - 203469.n_hn_l + 26154.4n_hn_l^2 - 817.246n_hn_l^3 + 6.65319n_hn_l^4 - 0.0308418n_hn_l^5).} \quad (A18)$$

[1] M. Gell-Mann, R. J. Oakes, and B. Renner, Phys. Rev. **175**, 2195 (1968).

[2] S. Aoki et al. (Flavour Lattice Averaging Group), Eur. Phys. J. C **80**, 113 (2020).

[3] See e.g., H. G. Dosch and S. Narison, Phys. Lett. B **417**, 173 (1998); M. Jamin, Phys. Lett. B **538**, 71 (2002).

[4] S. Narison, Phys. Lett. B **738**, 346 (2014).

[5] M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Nucl. Phys. B**147**, 385 (1979).

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

[7] S. P. Klevansky, Rev. Mod. Phys. **64**, 649 (1992); T. Hatsuda and T. Kunihiro, Phys. Rep. **247**, 221 (1994).

[8] C. D. Roberts and S. M. Schmidt, Prog. Part. Nucl. Phys. **45**, S1 (2000); R. Alkofer and L. von Smekal, Phys. Rep. **353**, 281 (2001).

[9] P. Maris, C. D. Roberts, and P. C. Tandy, Phys. Lett. B **420**, 267 (1998); P. Maris and C. D. Roberts, Phys. Rev. C **56**, 3369 (1997).

[10] K. Langfeld, H. Markum, R. Pullirsch, C. D. Roberts, and S. M. Schmidt, Phys. Rev. C **67**, 065206 (2003).

[11] D. R. Campagnari, E. Ebadati, H. Reinhardt, and P. Vastag, Phys. Rev. D **94**, 074027 (2016).

[12] See e.g., D. Becirevic and V. Lubicz, Phys. Lett. B **600**, 83 (2004); P. Hernandez, K. Jansen, L. Lellouch, and H. Wittig, J. High Energy Phys. 07 (2001) 018; A. Duncan, S. Pernice, and J. Yoo, Phys. Rev. D **65**, 094509 (2002).

[13] E. Marinari, G. Parisi, and C. Rebbi, Phys. Rev. Lett. **47**, 1795 (1981).

[14] L. Del Debbio, L. Giusti, M. Lüscher, R. Petronzio, and N. Tantalo, J. High Energy Phys. 02 (2006) 011; L. Giusti and M. Lüscher, J. High Energy Phys. 03 (2009) 013; H. Fukaya, S. Aoki, T. W. Chiu, S. Hashimoto, T. Kaneko, J. Noaki, T. Onogi, and N. Yamada, Phys. Rev. Lett. **104**, 122002 (2010); **105**, 159901(E) (2010).

[15] G. P. Engel, L. Giusti, S. Lottini, and R. Sommer, Phys. Rev. D **91**, 054505 (2015).

[16] T. Banks and A. Casher, Nucl. Phys. B**169**, 103 (1980).

[17] H. Leutwyler and A. V. Smilga, Phys. Rev. D **46**, 5607 (1992).

[18] A. V. Smilga and J. Stern, Phys. Lett. B **318**, 531 (1993); K. Zyablyuk, J. High Energy Phys. 06 (2000) 025.

[19] See e.g., for a recent review, M. Faber and R. Höllwieser, Prog. Part. Nucl. Phys. **97**, 312 (2017).

[20] J. Gasser and H. Leutwyler, Ann. Phys. (N.Y.) **158**, 142 (1984); Nucl. Phys. B**250**, 465 (1985).

[21] J. Bijnens and T. A. Lahde, Phys. Rev. D **71**, 094502 (2005); J. Bijnens and G. Ecker, Annu. Rev. Nucl. Part. Sci. **64**, 149 (2014).

[22] S. Descotes-Genon, L. Girlanda, and J. Stern, J. High Energy Phys. 01 (2000) 041; S. Descotes-Genon and J. Stern, Phys. Lett. B **488**, 274 (2000); S. Descotes-Genon, N. H. Fuchs, L. Girlanda, and J. Stern, Eur. Phys. J. C **34**, 201 (2004); V. Bernard, S. Descotes-Genon, and G. Toucas,
J. High Energy Phys. 01 (2011) 107; V. Bernard, S. Descotes-Genon, and G. Toucas, J. High Energy Phys. 06 (2012) 051; M. Kolesar and J. Novotny, Eur. Phys. J. C **78**, 264 (2018).

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

[24] J.-L. Kneur and A. Neveu, Phys. Rev. D **85**, 014005 (2012).

[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 M. B. Pinto, Phys. Rev. Lett. **116**, 031601 (2016); Phys. Rev. D **92**, 116008 (2015).

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

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

[30] P. A. Baikov, K. G. Chetyrkin, and J. H. Kühn, Phys. Rev. Lett. **118**, 082002 (2017); T. Luthe, A. Maier, P. Marquard, and Y. Schröder, J. High Energy Phys. 07 (2016) 127; F. Herzog, B. Ruijl, T. Ueda, J. A. M. Vermaseren, and A. Vogt, J. High Energy Phys. 02 (2017) 090; T. Luthe, A. Maier, P. Marquard, and Y. Schröder, J. High Energy Phys. 10 (2017) 166; K. G. Chetyrkin, G. Falcioni, F. Herzog, and J. A. M. Vermaseren, J. High Energy Phys. 10 (2017) 179; 12 (2017) 006(A).

[31] P. A. Baikov, K. G. Chetyrkin, and J. H. Kühn, J. High Energy Phys. 10 (2014) 076; T. Luthe, A. Maier, P. Marquard, and Y. Schröder, J. High Energy Phys. 01 (2017) 081.

[32] P. A. Baikov and K. G. Chetyrkin, Proc. Sci. RADCOR2017 (**2018**) 025.

[33] V. I. Yukalov, Theor. Math. Phys. **28**, 652 (1976); W. E. Caswell, Ann. Phys. (Berlin) **123**, 153 (1979); I. G. Halliday and P. Suranyi, Phys. Lett. **85B**, 421 (1979); P. M. Stevenson, Phys. Rev. D **23**, 2916 (1981); Nucl. Phys. **B203**, 472 (1982); J. Killinbeck, J. Phys. A **14**, 1005 (1981); R. P. Feynman and H. Kleinert, Phys. Rev. A **34**, 5080 (1986); A. Okopinska, Phys. Rev. D **35**, 1835 (1987); 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. **B359**, 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); S. Gandhi and M. B. Pinto, Phys. Rev. D **46**, 2570 (1992); H. Yamada, Z. Phys. C **59**, 67 (1993); K. G. Klimenko, Z. Phys. C **60**, 677 (1993); A. N. Sissakian, 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); Phys. Rev. D **60**, 085001 (1999).

[34] C. Arvanitis, F. Geniet, M. Iacomi, J.-L. Kneur, and A. Neveu, Int. J. Mod. Phys. A **12**, 3307 (1997).

[35] C. Arvanitis, F. Geniet, J. L. Kneur, and A. Neveu, Phys. Lett. B **390**, 385 (1997).

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

[37] 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); J. Zinn-Justin, arXiv:1001.0675.

[38] R. Guida, K. Konishi, and H. Suzuki, Ann. Phys. (N.Y.) **241**, 152 (1995); **249**, 109 (1996).

[39] H. Kleinert, Mod. Phys. Lett. B **17**, 1011 (2003); B. Hamprecht and H. Kleinert, Phys. Rev. D **68**, 065001 (2003); B. Kastening, Phys. Rev. A **68**, 061601 (2003); **69**, 043613 (2004).

[40] J.-L. Kneur, A. Neveu, and M. B. Pinto, Phys. Rev. A **69**, 053624 (2004).

[41] K. G. Chetyrkin and V. P. Spiridonov, Sov. J. Nucl. Phys. **47**, 522 (1988);

[42] K. G. Chetyrkin and J. H. Kühn, Nucl. Phys. **B432**, 337 (1994); K. G. Chetyrkin and A. Maier, J. High Energy Phys. 01 (2010) 092.

[43] K. G. Chetyrkin and A. Maier (private communication).

[44] J. A. M. Vermaseren, S. A. Larin, and T. van Ritbergen, Phys. Lett. B **405**, 327 (1997); K. G. Chetyrkin, Nucl. Phys. **B710**, 499 (2005); M. Czakon, Nucl. Phys. **B710**, 485 (2005).

[45] A. Maier and P. Marquard (private communication).

[46] K. G. Chetyrkin, J. H. Kuhn, and M. Steinhauser, Nucl. Phys. **B505**, 40 (1997).

[47] K. G. Chetyrkin, R. Harlander, and M. Steinhauser, Phys. Rev. D **58**, 014012 (1998); A. Maier, P. Maierhofer, and P. Marquard, Nucl. Phys. **B797**, 218 (2008).

[48] K. G. Chetyrkin, J. H. Kuhn, and C. Sturm, Eur. Phys. J. C **48**, 107 (2006).

[49] C. Sturm, J. High Energy Phys. 09 (2008) 075.

[50] K. Nakayama, H. Fukaya, and S. Hashimoto, Phys. Rev. D **98**, 014501 (2018).

[51] See the QCD review chapter by J. Huston, K. Rabbertz, and G. Zanderighi in M. Tanabashi et al. (Particle Data Group), Phys. Rev. D **98**, 030001 (2018).

[52] D. J. Gross and A. Neveu, Phys. Rev. D **10**, 3235 (1974).

[53] K. G. Chetyrkin, B. A. Kniehl, and M. Steinhauser, Phys. Rev. Lett. **79**, 2184 (1997); Nucl. Phys. **B510**, 61 (1998); K. G. Chetyrkin, J. H Kühn, and C. Sturm, Nucl. Phys. **B744**, 121 (2006).

[54] K. G. Chetyrkin, J. H. Kuhn, and M. Steinhauser, Comput. Phys. Commun. **133**, 43 (2000); F. Herren and M. Steinhauser, Comput. Phys. Commun. **224**, 333 (2018).

[55] A. Bazavov et al. (MILC Collaboration), Proc. Sci., CD09 (2009) 007 [arXiv:0910.2966].

[56] H. Fukaya, S. Aoki, T. W. Chiu, S. Hashimoto, T. Kaneko, J. Noaki, T. Onogi, and N. Yamada (JLQCD and TWQCD Collaborations), Phys. Rev. D **83**, 074501 (2011).

[57] S. Aoki et al. (PACS-CS Collaboration), Phys. Rev. D **79**, 034503 (2009).

[58] C. T. H. Davies, K. Hornbostel, J. Komijani, J. Koponen, G. P. Lepage, A. T. Lytle, and C. McNeile (HPQCD Collaboration), Phys. Rev. D **100**, 034506 (2019).

[59] P. Fritzsch, F. Knechtli, B. Leder, M. Marinkovic, S. Schaefer, R. Sommer, and F. Virotta, Nucl. Phys. **B865**, 397 (2012).

[60] B. Blossier, Ph. Boucaud, F. De soto, V. Morenas, M. Gravina, O. Pène, and J. Rodríguez-Quintero (ETM Collaboration), Phys. Rev. D **82**, 034510 (2010); K. Jansen, F. Karbstein, A. Nagy, and M. Wagner (ETM Collaboration), J. High Energy Phys. 01 (2012) 025; F. Karbstein, M. Wagner, and M. Weber, Phys. Rev. D **98**, 114506 (2018).

[61] F. Karbstein, A. Peters, and M. Wagner, J. High Energy Phys. 09 (2014) 114.
[62] K. Cichy, E. Garcia-Ramos, and K. Jansen, J. High Energy Phys. 10 (2013) 175; R. Baron et al. (ETM Collaboration), J. High Energy Phys. 08 (2010) 097; B. B. Brandt, A. Jüttner, and H. Wittig, J. High Energy Phys. 11 (2013) 034.

[63] S. Capitani, M. Lüscher, R. Sommer, and H. Wittig, Nucl. Phys. **B544**, 669 (1999); **B582**, 762 (2000); M. Gockeler, R. Horsley, A. C. Irving, D. Pleiter, P. E. L. Rakow, G. Schierholz, and H. Stuben, Phys. Rev. D **73**, 014513 (2006); N. Brambilla, X. Garcia i Tormo, J. Soto, and A. Vairo, Phys. Rev. Lett. **105**, 212001 (2010); **108**, 269903(E) (2012); M. Kitazawa, T. Iritani, M. Asakawa, T. Hatsuda, and H. Suzuki, Phys. Rev. D **94**, 114512 (2016); K. I. Ishikawa,

I. Kanamori, Y. Murakami, A. Nakamura, M. Okawa, and R. Ueno, J. High Energy Phys. 12 (2017) 067.

[64] T. W. Chiu and T. H. Hsieh, Nucl. Phys. **B673**, 217 (2003).

[65] C. T. H. Davies, C. McNeile, A. Bazavov, R. J. Dowdall, K. Hornbostel, G. P. Lepage, and H. Trottier, Proc. Sci. ConfinementX (2012) 042.

[66] C. A. Dominguez, N. F. Nasrallah, R. Rontsch, and K. Schilcher, J. High Energy Phys. 05 (2008) 020.

[67] S. Narison and R. Albuquerque, Phys. Lett. B **694**, 217 (2010); R. M. Albuquerque, S. Narison, and M. Nielsen, Phys. Lett. B **684**, 236 (2010).

[68] P. Gubler and D. Satow, Prog. Part. Nucl. Phys. **106**, 1 (2019).