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

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 {\em 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-loops RG properties. The RGOPT results overall show a very good stability and convergence, giving primarily the RG invariant condensate, $\langle \bar q q\rangle^{1/3}_{RGI}(n_f=0) = -(0.840_{-0.016}^{+0.020}) \bar\Lambda_0 $, $\langle\bar q q\rangle^{1/3}_{RGI}(n_f=2) = -(0.781_{-0.009}^{+0.019}) \bar\Lambda_2 $, $\langle\bar q q\rangle^{1/3}_{RGI}(n_f=3) = -(0.751_{-.010}^{+0.019}) \bar\Lambda_3 $, where $\bar\Lambda_{n_f}$ is the basic QCD scale in the \overline{MS} scheme for $n_f$ quark flavors, and the range spanned is our rather conservative estimated theoretical error. This leads {\it e.g.} to $ \langle\bar q q\rangle^{1/3}_{n_f=3}(2\, {\rm GeV}) = -(273^{+7}_{-4}\pm 13)$ 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.


I. INTRODUCTION
The chiral quark condensate qq is a main order parameter of spontaneous chiral symmetry breaking, SU (n f ) L × SU (n f ) R → 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: relates the condensate with the pion mass m π and decay constant F π 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 O(m 2 q ). Thus more direct "first principle" 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][9][10][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][14][15], directly related to the quark condensate via the Banks-Casher relation [16][17][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.
Our renormalization group optimized perturbation (RGOPT) approach [23][24][25] provides analytic sequences of (variational) nonperturbative approximations, having a non-trivial 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 the vacuum from the related spectral density, by going at the complete five-loop and partial six-loop level of our approximation.
In section 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 section 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 section 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 welldefined 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 arbitrary number of quark flavors can also be useful for other purposes irrespectively of our variational approach, most typically for perturbative matching of lattice simulation results. Section V give 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 Section VI we compare with other recent determinations, mainly from lattice simulations. Finally section VII presents a summary and conclusions, and an 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 ρ(λ) of the Dirac operator with the chiral condensate qq 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 (OPE) 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 qq (m) = −2 m ∞ 0 dλ ρ(λ) λ 2 + m 2 . (2.1) Eq.(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→0 qq = −πρ(0) , (2.2) if the spectral density at the origin can be determined. Note also that for non-zero fermion mass m, the spectral density is thus determined by the discontinuity of qq (m) across the imaginay axis: For nonvanishing quark mass m, qq has a nontrivial perturbative series expansion, ∼ m 3 f [ln(m 2 /µ 2 )], and its discontinuities are simply given by those coming from the perturbative logarithmic mass dependence. Therefore the above relation (2.3) also allows to calculate the corresponding perturbative spectral density. However, the λ → 0 limit, relevant for the true chiral condensate, trivially leads to a vanishing result, since perturbatively ρ(λ) ∼ λ 3 . But as we recall below a crucial feature of the variational RGOPT method is to circumvent this, giving a nontrivial result for λ → 0.

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 < δ < 1, interpolating between the (massive) free Lagrangian L f ree and the original (massless) Lagrangian L int respectively. This amounts first to the prescription: within some given (renormalized) perturbative expansion of a physical quantity P (m, g) (here g ≡ 4πα 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 δ at order k, the so-called δ-expansion [33], and afterwards δ → 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: determining a nontrivial dressed massm(g). The prescription is consistent with renormalizability [34][35][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 δ-expansion mainly for simplicity. In contrast we combine [23][24][25] the OPT Eq.(3.2) with renormalization group (RG) properties, by requiring the (δ-modified) expansion to satisfy, in addition to Eq.(3.2), a perturbative RG equation: where the (homogeneous) RG operator acting on a physical quantity is defined as 1 Note that once combined with Eq. (3.2), the RG equation takes a reduced massless form: 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 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]: Therefore Eqs. (3.5), (3.6) and (3.2) together completely fix optimized m ≡m and g ≡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 δ orders at least one of both the RG and OPT solutionsg(m) continuously matches the standard perturbative RG behavior for g → 0 (i.e. Asymptotic Freedom (AF) for QCD): 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 occurence of complex solutions is merely a consequence of solving exactly the (polynomial) Eqs.(3.2, (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 a 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 ρ(λ, a s ) was derived in [26] and is completely analogous to the standard RG equation, but with the mass replaced by the spectral parameter, One can next proceed to the modification of the resulting perturbative series ρ(λ, a s ) as implied by the δ-expansion, now, from Eq. (3.8) clearly applied not on the original mass but on the spectral value 3 λ: Consequently the mass optimization on qq thus translates into an optimization of the spectral density with respect to λ, at successive δ k order (see [26] for more details). Finally, as one last subtlety, note that the interpolation exponent a in Eq.(3.6) is universal in so far as the original expansion to be modified is itself (perturbatively) RG invariant. Now since m qq is the RG invariant quantity, rather than qq , 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 (3.11) which also maintains the occurence of essentially unique AF-matching solutions with a behavior similar to (3.7) (with m → λ understood).
2 A connection of the exponent a with RG anomalous dimensions/critical exponents had also been established previously in the D = 3 Φ 4 model for the Bose-Einstein condensate (BEC) critical temperature shift by two independent OPT approaches [39,40]. 3 We simplify notations with λ ≡ |λ| since it is necessarily positive. 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 qq up to six-loop order reads formally: where L m ≡ ln(m/µ), m ≡ m(µ) and a s ≡ α S (µ)/π in the M S scheme with renormalization scale µ. The two-loop contributions were calculated in the M S-scheme long ago, first in [41] (see also [36]). At higher k-loop orders (k ≥ 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 cancelled 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 qq 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 where the first term is the (homogeneous) RG operator given in Eq.(3.4) and Γ 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 Appendix). Therefore note that the RG consistency expressed by requiring Eq.(4.2) to hold perturbatively order by order, allows to determine all the logarithmic (L m ) p coefficients c ki , with k ≥ i + 2 at perturbative orders k from lower (< k) order coefficients and RG β(a s ) and γ m (a s ) functions up to order k − 2 and k − 1 respectively. In addition the knowledge of Γ 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 ≤ 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) B. Exact versus approximate determinations of c44 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+1 k single logarithm coefficient via RG properties, as explicited in Eq.(4.6) (see also Eq.(A12) in 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 where Π s (q 2 ) ≡ i d 4 xe iq.x 0|T J s (x)J s (0)|0 is the two-point scalar correlation function (the scalar current being defined as J s =qq). 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 Π 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 Π s (0) at three and four loops they are nonvanishing only for massive quark contributions ∝ n h ). The nonsinglet four-loop nonlogarithmic contribution to Π s (0) is given in the M S-scheme in Eq.(B.1) of [49], a result that we recast here for completeness in our normalization conventions: 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 Π s (0). Accordingly from a straightforward integration from Eq.(4.7) with input (4.8), we obtain the 'incompletesinglet' (IS) approximation of c 44 6 : Note that the first (dominant) term in Eq.(4.10) is the pure gauge contribution, while terms ∝ 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. 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 | qq | 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], where we define with coefficients determined order by order by 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 β(a s ) and γ 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 Appendix A (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 pertubative 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 non-logarithmic terms to zero, those having obviously no discontinuities, while replacing all powers of logarithms, using m → i|λ| etc., as leading to the following substitution rules for the first few terms where L λ ≡ ln(λ/µ) (note the π 2k terms appearing starting at order ln 3 m).
We obtain in this way the perturbative spectral density up to six-loop order formally: where the coefficients ρ ki for k ≥ 3 are straightforwardly related to the c ki of the original condensate using (4.16) as follows: and so on at higher (six-loop) order (see Eq. (A15) in Appendix). At this stage, before proceeding with RGOPT, we remark that the above (ordinary) perturbative spectral density ρ(λ) expression for arbitrary n l , n h in Eq. (4.17) can be useful for different purposes, independently of the RGOPT approach. For instance it should allow 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 α 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, Eq.(3.10), Eq.(3.8), applied to the spectral density Eq.(4.17) at successive orders, after the modifications implied by Eq.(3.9), (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 selfcontainedness, 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 λ contribution, one obtains for Eq (3.10) the trivial optimized solution, λ = 0. Thus nontrivial solutions occur starting at next-to-leading (NLO) two-loop order of the modified perturbation. Accordingly at NLO order the modifed series reads  Tables I, II for n f = 2, 3 respectively, using also (2.2). These results used the RG Eq. (3.8) at one-loop order, that give 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 effects. Doing this gives the results in the second lines of Tables I, II for n f = 2, 3. Those results, to be considered more accurate, show a substantial decrease of the optimal coupling α 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Λ 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,Λ = µe −1/(2b0g) , consistently 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Λ (see, eg., [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Λ. The four-loop QCD scaleΛ expression reads in our normalizations (where g ≡ 4πα S (µ)): 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 qq RGI , which can be more appropriately compared between different perturbative orders. It is defined in our normalization as where higher-order terms not shown here are easily derived, since only depending on the RG coefficients b i , γ i known up to five loops. Remark, however, that our RGOPT optimization also fixes a scale, simply obtained from using Eq.(5.2) (or its lower order equivalent) forΛ(g), that are given indicatively in Table I, II. We stress that the optimal couplingα S and corresponding optimal scaleμ, or the optimal spectral parameterλ, 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α S andλ within the quantity being optimized, here ρ(λ, α 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 2 s , order, the n f dependence appears explicitly within the perturbative coefficients of the spectral density, see Fig 1 and the last a 2 s coefficient in Eq. (4.17). There occurs two real solutions forL λ ,α S , but the selection of the unique physical solution is unambiguous since only one is clearly compatible with AF behavior for (1), both for the RG and OPT equations. In contrast the other real solution has for g → 0 a coefficient of opposite sign to AF, and gives lnλ/µ > 0, which is incompatible with perturbativity, since we expect µ ≫λ similarly to the perturbative range µ ≫m ∼Λ for the original expansion with mass dependence. As stressed above in Sec. III the occurence 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, II respectively 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λ,α S but there are no possible ambiguities since once more all solutions are eliminated from the AF-matching requirement, except a single one, with α S > 0 andL λ < 0 as expected.
One observes a further decrease of the optimal couplingα S from three to four loops to more perturbative values, as well as the corresponding decrease ofL λ , meaning thatμ is also larger. The stabilization/convergence of the results is clear for the scale-invariant condensate qq RGI given in Tables I, II, which at four-loop order has almost no variation upon RG equation truncations 9 . Note that the optimal valuesα S decreases substantially with increasing orders as compared to the lowest nontrivial order result above, thus indicating more perturbatively reliable results, moreoverα S appears to somehow stabilize at three and four loops. Notice also that compared with the more than 10% change inα S upon going from two to three loops, the final physical condensate value only varies by 0.25%, showing a strong stability. Also, while qq 1/3 /Λ changes by about 20% compared to the crude two-loop result in the first lines of Tables, 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 four loops, respectively as L RG λ (α S , k), L OP T λ (α 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 λ values for any α S > 0 in Fig. 2) are real for any α S > 0 and intersect at a unique value, that determines unambiguously the solution, compare 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 unconveniently become complex (conjugate)-valued within the range 0.25 α S 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 β and γ m functions within the optimization, after performing consistently the δ-expansion in Eq. (3.9) to order δ 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 δ 5 ). More precisely all the presently known five-loop RG coefficients, together with the complete four-loop coefficients, allow to determine exactly the six-loop coefficients c 6k , 0 ≤ k ≤ 4 of ln 6−k (m/µ) 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 Appendix. Consequently from Eq.(4.16) all the six-loop logarithmic terms, ρ 6k L 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 occurence of rather unwelcome instabilities for the strictly five-loop results.

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λ = −0.5699,α S = 0.5963, which gives qq  Table I. A directly related issue is the anomalously large optimized coupling that corresponds to this solution, α S ∼ 0.6, in contrast with the regularly decreasing coupling obtained at increasing orders up to four loops, in Tables I, 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 Eqs. respectively, that both give quartic equations in L λ ≡ ln(λ/µ) at five-loops (thus respectively giving L RG λ (α S , k), L OP T λ (α S , k) with k = 1, ..4), and examining which one(s) exhibit the perturbative AF-matching behavior, and whether the latter are matching the optimized Lλ,α 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 λ and α S ): in contrast with the four-loop results in Fig. 2, the RG branches now also become complex, similarly to the OPT branches, within a rather important α S range: 0.27 α S 0.56, and the only real intersection occuring at α S ≃ 0.596, Lλ ≃ −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 Γ 0 4 in Eq.(A5) with Γ 0 3 in Eq.(A4) in Appendix). (In contrast, as illustrated below, going from the four-loop to the five-loop β-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 β(α S ), γ m (α S ) functions). Note that Γ 0 4 enters the five-loop condensate L m coefficient c 54 (n f ), but it is not the sole contribution: the net effect from Γ 0 4 is typically that |c 54 (n f )/c 50 (n f )| ∼ 10 roughly, which may be compared qualitatively with the similar four-loop quantities giving |c 43 /c 40 | ∼ 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 ∼ π 4 , appearing for the first time at five-loop order in ρ(λ, a s ), within the nonlogarithmic coefficient. More precisely, it is the last term of Eq.(4.19), 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 on Fig. 3 that the AF-matching RG and OPT branches have not disappeared but just became complex (conjugate) valued within a certain α 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λ ≃ −0.778 ± 0.303i,α S ≃ 0.358 ± 0.0537i, that gives for the (RG invariant) condensate: ≃ −(0.801 ± 0.0195i)Λ(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α 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 occurence of nonreal solutions is only a consequence of solving exactly the RG and OPT polynomial equations in L λ , 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 SubSec. V B). Let us examine in turn those two possibilities.

Perturbatively truncated five-loop RG solutions
At five loops, instead of solving exactly the relevant RG and/or OPT optimization Eqs.(3.8), (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 (α 4 S ) spectral density Eq.(4.17) involves up to α 9 S terms, due to the highest five-loop RG contributions ∝ b 4 α 6 S , γ 4 α 5 S respectively. But b 4 , γ 4 appear first at order α 5 S , α 4 S respectively. Accordingly a presumably sensible procedure is to truncate [24] the RG equation, suppressing higher-order terms in α 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 , γ 5 (presently not known), would enter first respectively the α 6 S , α 5 S coefficients of the RG equation. Therefore it appears sensible to truncate any α k S , k ≥ 6 in the result of Eq. (3.8), that would be anyway affected by presently unknown higher orders. Further truncating the α 5 S term implies, however, losing any dependence from the five-loop b 4 (while it still involves the five-loop γ 4 one). Accordingly we found instructive to consider the effects of successive truncations, progressively suppressing the highest α 9 S down to α 6 S (or even possibly α 5 S ) terms and comparing. This is done below, with all results compiled in Table III obtained by optimizing the spectral density ρ(λ, α S ) and keeping only the AF-matching branch solution (unique at a given order).
TABLE III. n f = 2 results at RGOPT δ 4 (five-loops) and partial δ 5 (six-loops) for the (RG invariant) condensate qq 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 Im[Lλ] 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α S value are also more reasonably perturbative, and close to the real four-loop results of Table I. Truncating maximally the RG equation (namely by all terms α k≥6 S , 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 qq 1/3 RGI ≃ −(0.800 ± 0.0005i)Λ 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.

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 (λ) coefficients given explicitly in Appendix (see Eqs.(A14), (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 ∝ π 4 (see the last terms in Eq.(4.16)), that tend to partially balance the instability triggered by π 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 towards the exact result (known for the GN model), those solutions being destabilized each time novel π 2k contributions appear first, at increasing orders 10 . However, if keeping only fixed π 2k terms (namely discarding π 2k+2 , etc, terms appearing at higher orders), remarkably at sufficiently high fixed order all the π 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 ρ(λ) (all except the nonlogarithmic one, ρ 66 in Eq. (4.17). Accordingly one should keep in mind that it remains an approximation to the complete six-loop results, since ρ 66 involves the presently unknown six-loop vacuum anomalous dimension and the five-loop nonlogarithmic coefficient c 55 . Therefore we simply set ρ 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 α S range, and their intersection occur for a substantially smallerα 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 λ , α 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 π 4 terms occuring first at five loops.

Summary of n f = 2 results
As a tentative summary of the previous n f = 2 investigation: • At five loops, the impact of both large five-loop vacuum energy anomalous dimensions and (more importantly) the first occurence of π 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 unmandatory 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 α k≥6 S , see Table III) and can be included within the theoretical uncertainties. • Next, going to six loops restores a real unique AF-matching solution, that results from a partial balance of the destabilizing π 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 similarly 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λ ≃ −0.693411,α S ≃ 0.598. This gives for the RG invariant condensate: ≃ −0.813Λ 3 , thus a large shift from four-loop results of Table II. But again upon examining the AF branches these do not match with this real solution, and the correct AF-matching but complex-valued branch gives a more reasonable result, with small imaginary parts and real part closer to four-loop real results (see Table IV).
Similarly 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 much similar to n f = 2: at five loops one fails to recover strictly real AF-matching solutions, but the maximal truncation (discarding α k≥6 S , not loosing 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, similarly 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 subsection IV B and relevant for n f = 2, 3, of using for the fourloop nonlogarithmic 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 subSec.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 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 similarly 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 from using the exact c 44 in Tables III, 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, IV by about ∼ 1% (∼ 2%) lower inΛ units for n f = 2 (n f = 3) respectively. All these features are consistent with the fact that the loss of real solution at five loops is essentially due to the occurence of relatively large π 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 | qq | 1/3 that is about ∼ 1%(2%) smaller in magnitude respectively for n f = 2 (n f = 3), which again reflects a good overall stability.  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 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 | qq | 1/3 is driven to about ∼ 2.5% higher values when going from four to five loops, as similarly observed above for n f = 2, 3 in Tables III, IV. As above mentioned this is essentially traced to the instability from the first occurence of π 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  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 crosscheck of our results. We have explored variants similarly to other n f values above but simply summarize here the main results. At four-loop order one obtains the unique real solution: Similarly to n f = 0, 2, 3 cases, once more the five-loop results produce a substantial ∼ 2.8% increase of the condensate as compared to 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 five-loop results. Note indeed that if discarding the latter, the combined 4-loop and 6-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 qq 1/3 where we give only a three digits accuracy given the uncertainties. Eqs (5.7),(5.8), (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 µ ≃ 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 α s (and other related) uncertainties.

VI. qq (µ = 2 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 µ ′ , using again (5.3) at five-loop order, now taking g ≡ 4πα S (µ ′ ), after evolving α S (µ) at five-loop order of Eq.(5.2) towards the conventional scale µ ′ = 2 GeV. The overall reliability of this (perturbative) evolution is to be assessed on the ground that the primary RGOPT results above at four-, five-and six-loops are obtained at reasonably perturbative optimized scale values (3.3Λ MS μ 4.8Λ 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 α S determinations in the perturbative range. We also account properly for the charm quark mass threshold effects [53] on α S (µ ∼ m c ). From the most recent world average α S (m Z ) value [51]: where the first error is our rather conservative theoretical RGOPT uncertainty from Eq.(6.4) and the second one is fromΛ(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 α S world average with substantial increase of uncertainties (see [51] for detailed explanations on these features).
To compare with other independent determinations, first the most precise n f = 3 lattice determination we are aware of, in the chiral limit, is qq [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 ± 6 ± 24 [56] to 290 ± 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 heavy-strange pseudoscalar correlator, has obtained ss 1/3 (2 GeV) ∼ −(296 ± 11) 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Λ(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 α S as above. Nevertheless, given that the (optimized) coupling values obtained in Table I, 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Λ(n f = 2). To our knowledge there are not so many nonperturbative results for Λ(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:Λ (n f = 2) = (310 ± 20) MeV. (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] 12 . After RG evolution up to 2 GeV we obtain accordingly from Eq. where again the first error range is our RGOPT uncertainty from Eq.(6.7) while the second one is from theΛ (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]: qq 1/3 n f =2 (µ = 2GeV) = −(261 ± 6 ± 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]: | qq | 1/3 n f =2 = −(266 ± 10) MeV, where the uncertainties include both systematic and statistical ones.
One may also compare with recent results from spectral sum rules [4]: ūu 1/3 (≡ d d 1/3 )(2 GeV) ≃ −(276 ± 7) MeV. But keeping in mind that the latter sum rules actually determine precisely the current quark masses, so that the ūu 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Λ (2) is not yet at the level ofΛ(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Λ(2): taking the above quoted estimate of the condensate given by [2], this gives Λ(n f = 2) = (308 +4 −6 ± 12)MeV. (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 α S (µ) and use Eq. (5.3) at five loops but in the appropriate n f = 0 approximation. As previously we needΛ(n f = 0), available from various different approaches with lattice simulations. We rely on the average performed in [2], combining different precise lattice results [63]:Λ (n f = 0) = (257 ± 7) MeV. (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.472fm, which for n f = 0 amounts merely to a defining convention forΛ(n f = 0). Next the RG evolution from Eq. (6.12) leads to 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 crosscheck of our method, since the quenched approximation is anyway not very realistic.
D. Further discussion on qq (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Λ 3 has a sizable but moderate dependence on the number of flavors n f : there is a clear trend that | qq 1/3 |(n f )/Λ(n f ) decreases regularly, roughly linearly by about 4% for n f → n f + 1 (that is clear at least from the studied n f ≤ 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Λ(n f ), that appears rather to increase with n f for n f ≤ 3 (at least if considering the low latticeΛ(0) value (6.10, but it is not so clear from comparingΛ (2) andΛ (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]: ss / ūu = 0.74 +0. 34 −0.12 , (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 M S 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), (6.5) are further affected by present uncertainties in perturbative evolution andΛ values. Eqs. (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 occurence at growing orders of new large π 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 | qq | 1/3 RGI (n f ) obtained at five-loop RGOPT order are largely an artifact of the instability from the π 4 discontinuity terms appearing first at five-loop order. Nevertheless we incorporate more conservatively the differences between four-, five-and sixloop 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), (6.8) are more affected by the present uncertainties on the basic QCD scaleΛ, both for n f = 2 and n f = 3. (To possibly get rid ofΛ uncertainties, particularly for n f = 2, one could in principle apply RGOPT directly to a more physical RG invariant quantity, like qq RGI /F 3 π , combining the present analysis with the one in [25] for F π : 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 f values. Our results for n f = 3 are perhaps of particular interest, given that other independent determinations are either not in the chiral limit or, concerning lattice results, are affected by still rather important uncertainties in the chiral extrapolation [2]. Finally our results indicate a moderate flavor dependence of the qq 1/3 n f /Λ n f values in the chiral limit for n f ≤ 3.