Calculation of the non-perturbative strong coupling from first principles

The success of the screened massive expansion is investigated in the framework of a screened momentum-subtraction scheme for the running of the strong coupling in pure Yang-Mills theory. By the exact Slavnov-Taylor and Nielsen identities, a very predictive and self-contained set of stationary conditions are derived for the optimization of the fixed-coupling expansion, yielding explicit analytical one-loop expressions for the propagators, the coupling and the beta function, from first principles. An excellent agreement is found with the lattice data. In the proposed screened renormalization scheme, a monotonic running coupling emerges which saturates in the IR at the finite IR stable fixed point $g=9.40$ where the beta function crosses the zero. A simple analytical expression is derived for the leading behavior of the beta in the IR.

In the last years, a purely analytical approach to the exact gauge-fixed Lagrangian of QCD has been developed by a screened massive expansion [47][48][49][50][51][52]. Preliminary variational calculations [31][32][33] had shown that PT might actually work at low energy and give reasonable results at the lowest orders of approximation if the expansion point is changed and the expansion is taken around a massive gluon propagator. Regardless of the accuracy of the zeroth order gluon propagator, provided that it has a mass scale, the higher-order terms become small, suggesting that the failure of PT might be a consequence of the bad choice of expanding around the usual massless gluon propagator of the high energy theory. On the other hand, the remarkable discovery of Tissier and Wschebor [44,45], that PT is viable if a gluon mass is added by hand to the Lagrangian, suggested that the best expansion point might actually be the simple free massive gluon which emerges by just adding a gluon mass term to the quadratic part of the Lagrangian and subtracting it again from the interaction, thus leaving the total action unchanged. Expanding around the modified massive quadratic part of the action, leads to the screened massive expansion which was developed in Refs. [47,48] and has the merit of providing explicit and very accurate analytical expressions for the propagators of the gauge-fixed Faddeev-Popov Lagrangian, without adding any phenomenological parameter [49]. The method was extended to the full QCD in Ref. [50], to a generic covariant gauge in Ref. [51] and to finite temperature in Refs. [52,53], enforcing the idea that most of the nonperturbative effects can be embedded in the gluon-mass parameter [54][55][56][57].
A great advantage of the screened expansion is that, when optimized, it provides explicit analytical expressions which can be continued to the whole complex plane, gaining access to information that cannot be currently obtained by other theoretical tools [58][59][60][61]. At variance with the Curci-Ferrari model which was studied in Refs. [44,45] as a low energy effective model, in the screened expansion no mass term is added to the total action which is not modified at all and still has the usual exact Becchi-Rouet-Stora-Tyutin (BRST) symmetry. It is the modified expansion point which breaks the symmetry softly at any finite order, so that the exact constraints inherited by the BRST symmetry, like Slavnov-Taylor (ST) and Nielsen identities [62], can be used as a self-contained criterion for optimizing the expansion and fixing the finite part of the renormalization constants, which are usually scheme-dependent and undetermined in PT. Without any phenomenological parameter left, the optimized expansion is highly predictive, as shown in Ref. [51] where the optimal gluon propagator was determined by the Nielsen identities.
In this paper, the optimization of the screened expansion for pure Yang-Mills theory is revised and extended to the ghost propagator, which plays a key role for determining the running coupling in the Landau gauge [63,64]. A natural variant of the momentum-subtraction (MOM) scheme is proposed by requiring that, at the subtraction point, the renormalized gluon propagator is equal to the screened massive free-propagator instead of the massless propagator of standard PT. In that screened scheme, the expansion is optimized using Nielsen identities, ST identities and enforcing the multiplicative renormalization of the theory. From first principles, without any phenomenological parameter, the gluon and ghost propagators are found in excellent agreement with the lattice data. Moreover an ab initio lattice-independent measure of the strong coupling is provided in the continuum, avoiding the effects of a finite lattice spacing [18].
In the screened MOM scheme, the strong coupling is a monotonic function which saturates in the IR at a finite fixed point where the beta function meets a zero. Thus, the proposed scheme could have more phenomenological relevance than the standard MOM scheme. Furthermore, analytical expressions can be easily derived for the strong coupling and the beta as functions of the renormalization scale. The non-perturbative beta function is fully characterized and a simple explicit expression is provided for the leading behavior at the finite IR fixed point. The paper is organized as follows: in Section II, the renormalization of Yang-Mills theory is briefly reviewed in order to fix the notation and recall some general results; in Section III a screened MOM scheme is introduced for the screened massive expansion; in Section IV the fixed-coupling expansion is optimized by a set of stationary conditions which derive from the BRST symmetry of the total action; in Section V the predictions of the optimized expansion are compared with the available lattice data and the properties of the beta function are described in detail; finally, a brief discussion of the main results is given in Section VI. Explicit analytical expressions and the leading behavior of the propagators, strong coupling and beta function are collected in Appendix A.

II. RENORMALIZATION OF YANG-MILLS THEORY
In a linear covariant ξ-gauge, the exact gauge-fixed bare Lagrangian of pure Yang-Mills SU(N) theory can be written as and L F P is the ghost term arising from the Faddeev-Popov (FP) determinant. Here, the tensor operator iŝ and the gauge field operators satisfy the SU (N ) algebrâ The total action can be splitted as S tot = S 0 +S I where the quadratic part is while the interaction contains the three vertices In Eq. (5), ∆ 0 and G 0 are the standard free-particle propagators for gluons and ghosts and their Fourier transforms are having used the transverse and longitudinal projectors In the present paper, we will take ξ → 0 and use the Landau gauge which is a Renormalization Group (RG) fixed point and is the most studied gauge on the lattice. Moreover, we will take advantage of Taylor's nonrenormalization theorem [63] which is only valid in the Landau gauge. Therefore, we will not bother about the gauge parameter ξ in the following discussion.
In the above equations, the fields and the coupling must be regarded as the bare A B , c B , g B . The total Lagrangian in Eq.(1) has a BRST symmetry which can be used for proving exact Slavnov-Taylor (ST) identities, Nielsen identities, and the multiplicative renormalizability of the theory. Then, we can follow the standard path of renormalization and introduce usual renormalized fields and coupling A, c, g having dropped obvious color indices. In terms of the renormalized fields, the Lagrangian can be written in the same identical form as in Eq.(1), provided that a set of counterterms are added to the interactions of Eq. (7). For instance, with an obvious shorthand notation, the gluon quadratic part becomes , so that the second term on the right hand side can be regarded as a counterterm which contributes to the polarization at tree-level.
On the same footing, the gluon-ghost vertex L ccg can be written as where the vertex renormalization constant Z c 1 is defined as and δZ c 1 = (Z c 1 − 1) is the prefactor of the vertex counterterm in the last term of Eq. (12).
An important consequence of ST identities is that all divergences can be eliminated by the same Z g for all vertices, since the ghost-gluon, three-gluon and four-gluon vertex renormalization constants, Z c 1 , Z 3g 1 and Z 4g 1 respectively, must satisfy the exact equations In perturbation theory (PT), the diverging parts of the renormalization constants satisfy Eqs. (14) at each order, but their finite parts might not respect the same constraints. Actually, in PT the finite parts of the renormalization constants are not determined at all and depend on the renormalization scheme. While they are simply ignored in the Minimal Subtraction scheme (M S), a more physical criterion is provided by the Momentum Subtraction scheme (MOM). The optimization of that choice leads to a better convergence of PT if the finite parts can be used as variational parameters, yielding an optimized PT by variation of the renormalization scheme [65]. Of course, the use of that method requires the existence of a function which is known to be minimal when the approximate result approaches the exact one, or at least stationary, according to Stevenson's principle of minimal sensitivity [66].
As discussed in Ref. [51], one important merit of the screened massive expansion arises from the apparent drawback of having splitted the action in two parts that are not BRST invariant. The exact constraints arising from BRST invariance cannot be satisfied exactly at any finite order of the expansion, but they must be satisfied by the exact result. Thus, ST and Nielsen identities provide a criterion for optimizing the expansion and fixing the finite part of the renormalization constants. In Ref. [51] the optimal gluon propagator was determined by requiring that the poles and the residues are gauge parameter independent, as required by Nielsen identities. At one-loop, the method provides a very accurate analytical expression for the gluon propagator, in very good agreement with the lattice data.
The method can be extended to the ghost propagator, which plays a key role for determining the running coupling because of the non-renormalization of the ghostgluon vertex in the Landau gauge. The extension must rely on ST identities and on the multiplicative renormalization of the theory, according to Eqs. (14).

III. SCREENED MOMENTUM-SUBTRACTION
A screened expansion for the exact renormalized Yang-Mills (YM) Lagrangian was developed in Refs. [47,48] and extended later to finite temperature in Refs. [52,53] and to the full QCD in Ref. [50], where a set of chiral quarks were included. The extension to a generic covariant gauge [51] has already shown the predictive power of the method when the expansion is optimized by the constraints of BRST symmetry.
The screened expansion arises by changing the expansion point after having renormalized the fields and the coupling as discussed above. Thus, the exact renormalization constants still satisfy Eqs. (14) because of the ST identities.
Following Refs. [48,51] the new screened expansion can be recovered by just adding a transverse mass term to the quadratic part of the action and subtracting it again from the interaction, leaving the total action unchanged. After having renormalized the fields and the coupling, we add and subtract the action term where the vertex function δΓ is a shift of the inverse propagator and ∆ m µν is a new massive free-particle propagator Adding that term is equivalent to substituting the new massive propagator ∆ m µν for the old massless one ∆ 0 µν in the quadratic part. Thus, the new expansion point is a massive free-particle propagator for the gluon, which is much closer to the exact propagator in the IR. In order to leave the total action unaffected by the change, the same term is added in the interaction, providing a new interaction vertex δΓ, a two-point vertex which can be regarded as a new counterterm. Dropping all color indices in the diagonal matrices and inserting Eq. (8) and (17) in Eq.(16) the vertex is just the transverse mass shift of the quadratic part δΓ µν (p) = m 2 t µν (p) (18) and must be added to the standard set of vertices arising from Eq.(7). The new vertex does not contain any renormalization constant and is part of the interaction even if it does not depend on the coupling. The expansion must be regarded as a δ-expansion more than a loop expansion, since different powers of the coupling coexist at each order in powers of the total interaction. Then, the proper gluon polarization and ghost self energy can be evaluated, order by order, by PT. In all Feynman graphs, the internal gluon lines are replaced by the massive free-particle propagator ∆ m µν and new insertions must be considered of the (transverse) two-point vertex δΓ µν . For further details of the screened expansion we refer to Refs. [48,51] where explicit analytical expressions are reported up to third order in the δ-expansion and one-loop.
Since the total gauge-fixed FP Lagrangian is not modified and because of gauge invariance, the longitudinal polarization is known exactly and is zero. We can write the exact polarization as so that, in the Landau gauge, the exact gluon propagator is transverse and defined by the single scalar function ∆(p).
From now on we switch to the Euclidean formalism and denote by p 2 the Euclidean square which is negative on the positive real axis of Minkowski space. The exact (dressed) gluon and ghost propagators are given by the functions where the proper gluon polarization Π and ghost selfenergy Σ are the sum of all one-particle-irreducible (1PI) graphs in the screened expansion, including all counterterms.
We observe that the mass parameter m is arbitrary and is the only energy scale in the expansion. Since there are no other energy scales in the Lagrangian, we can take m as the natural choice of energy units. Like for Λ QCD or for the inverse of the spacing on a lattice, the actual value of m can only be determined by the phenomenology. We will fix the scale m by comparison with the lattice data of Ref. [18].
Of course, any renormalization of m does not make any sense, since it is an arbitrary quantity which is added and subtracted again in the action. Rather, like Λ QCD , it must be regarded as a RG invariant which fixes the physical energy scale. On the other hand, BRST invariance protects the expansion from the appearing of any spurious diverging mass term, so that no mass counterterm is required and no diverging renormalization constant would require a renormalization of m. Actually, the added mass term breaks the BRST symmetry of the quadratic part S 0 and of the interaction S I when taken apart. Therefore, many constraints arising from BRST are not satisfied exactly at any finite order of the screened expansion. While the soft breaking has no effect on the UV behavior and on the diverging parts of the renormalization constants, some spurious diverging mass terms do appear in the expansion at some stage. However, as discussed in Refs. [48,50,51], the insertions of the new vertex δΓ, Eq.(18), cancel the spurious divergences exactly, without the need of any mass renormalization, as a consequence of the unbroken BRST symmetry of the whole action S 0 + S I . That aspect makes the screened expansion very different from effective models like the Curci-Ferrari model where a bare mass term is present in the Lagrangian from the beginning.
The exact self energies can be written as where the tree-level contribution m 2 comes from the new two-point vertex δΓ in Eq.(18), while the tree-level terms −p 2 δZ A , p 2 δZ c arise from the respective field-strength renormalization counterterm, which was explicitly shown in Eq.(11) for the gluon. The proper functions Π loop , Σ loop are given by the sum of all 1PI graphs containing loops.
The diverging parts of δZ A , δZ c cancel the UV divergences of Π loop and Σ loop , respectively. Since these divergences do not depend on masses, they are exactly the same as in the standard PT, so that in the M S scheme Z A and Z c have their standard expressions, as manifest in the explicit one-loop calculation [48]. The finite parts of δZ A , δZ c are arbitrary and depend on the scheme. On the other hand, the self energies contain an arbitrary term Cp 2 , where C is a constant which depends on the regularization method.
The MOM scheme is usually regarded as the most physical way to fix the finite parts for QCD in the Landau gauge. However, for the screened expansion, the most appropriate scheme seems to be a screened version of the standard MOM, which we call screened-MOM (SMOM). As suggested by Eq.(21), the SMOM scheme is defined by requiring that the dressed propagators are equal to the free-propagators of the screened expansion at the scale p = µ: at variance with the usual MOM scheme where the gluon propagator is set as ∆(µ) = 1/µ 2 . We observe that the validity of Taylor's non-renormalization theorem [63] is not affected by the masses and still holds in the screened expansion and in the Curci-Ferrari model [67,68]. Thus, we can set Z c 1 = 1 and require that the finite part of Z g satisfies the exact Eq.(13) which completes the definition of the scheme In Ref. [45], Eq.(23) was used in a MOM scheme for the Curci-Ferrari model, but Eq.(24) was different because the mass was a running coupling while here m is a RG invariant.
The SMOM scheme is non-perturbative and Eqs. (23), (24) are exact equations governing the multiplicative renormalization of the exact dressed propagators in Eq. (21). The approximate propagators that arise from the screened expansion might not be fully compatible with Eq. (24) which is expected to be satisfied in a perturbative sense, order by order. That is obviously true for the diverging parts which are the same as in the M S scheme and then satisfy Eq.(24) at each order of the screened expansion.
By a dimensional argument, the finite parts of the selfenergies can be written as power expansions where s = p 2 /m 2 and we use the short hand notation The functions F (s), G(s), F 2 (s), etc., are adimensional functions of the variable s and do not contain any parameter. Explicit analytical expressions for the one-loop functions F (s) and G(s) were evaluated in Ref. [48] up to third order in the δ-expansion and are reported in Appendix A. At one-loop, in the SMOM scheme, the subtracted self energies in Eq. (22) can be written as where t = µ 2 /m 2 is the renormalization scale in units of m. In the above equation, the normalization conditions of Eq.(23) are evident since Π(µ) = Σ(µ) = 0. Moreover, any arbitrary additive constant in the definition of the functions F (s), G(s) is subtracted and made irrelevant. By comparison with Eq. (22), the finite parts of the renormalization constants follow We observe that the two terms in δZ A are of different orders. The second one is of first order in the δ-expansion (i.e counting the number of vertices) but does not vanish in the limit α → 0 since it arises from the two-point vertex δΓ µν which adds a tree-level term m 2 to the gluon polarization in Eq. (22). Inserting in Eq.(21) and using a screened definition of the gluon dressing function J m (s) we can write the renormalized dressing functions at the scale t as where α(t) is the renormalized couplig α at the scale t.
It can be easily checked that J m (t, t) = χ m (t, t) = 1, as required by the SMOM scheme.
In the standard MOM scheme, the gluon dressing function is defined as J 0 (s) = p 2 ∆(p) and by comparison with the definition of J m (s), Eq. (29), the two functions are related by Multiplicative renormalization requires that we can define renormalization constants such that the ratios do not depend on the momentum s. Here, the last terms of the equalities arise by setting s = t ′ . Quite generally, the one-loop dressing functions can be multiplicatively renormalized, but in a very limited momentum range and for a very small coupling α ≪ 1. For instance, for the ghost case, in the limited range s ≈ t ≈ t ′ ≫ 1 we can set α = α(t) = α(t ′ ) + O(α 2 ) and using Eq.(30) That limited range can be extended numerically by evaluating RG improved dressing functions. Because of the special nature of the screened expansion, especially for the gluon dressing function, it is instructive to recover the RG integration by a direct calculation rather than using the anomalous dimensions. We will check that the same result is obtained by the anomalous dimensions if they are evaluated from the explicit one-loop renormalization constants of Eq. (28).
Using their definition, Eq.(30), the one-loop gluon dressing functions satisfy where δt = t ′ − t. Then, in the limit δt i = (t i+1 − t i ) → 0, setting t 1 = t and t N +1 = s, where F ′ (x) is the derivative of the function F (x). The second term can be integrated exactly yielding For the ghost dressing function, from Eq.(33) the same calculation gives We can easily check that, since dµ/µ = dt/(2t), Eq. (35) and Eq.(37) assume the more familiar shape where the anomalous dimensions γ A , γ c are given by the same explicit expressions which follow, by a direct calculation, from the renormalization constants of Eq.(28) We notice that, at one loop, the tree-level term 1/t must be retained in the inverse of (1 + δZ A ) since it is not of order O(α), yielding the same identical result found in Eq.(35) by a direct integration. The running coupling α(t), which is required for the numerical integration in Eqs. (36), (37), follows from ST identities which in the SMOM scheme, through Eqs. (24),(32) yield [64] which completes the definition of the RG improved oneloop dressing functions. The closed set of non-linear coupled integral equations, Eqs. (36), (37) and (40) can be solved numerically, starting from some initial point α(t 0 ). We are not pursuing the approach further since we are mainly interested in an analytical description of the IR limit, say below 2 GeV, where the fixed-coupling optimized screened expansion already provides excellent results. Actually, the main advantage of the fixed-coupling expansion is that it provides simple analytical expressions which can be continued to the whole complex plane in order to explore the analytic properties. A feature which makes the expansion a unique tool for accessing information that are not currently reachable by lattice calculations or other numerical methods. Eventually, the RG improved dressing functions might be useful for a matching between the optimized expansion and ordinary PT above 2 GeV, since the anomalous dimensions γ A , γ c in Eq. (39) tend to the standard result of ordinary PT in the limit t ≫ 1 (µ ≫ m).

IV. OPTIMIZED SCREENED EXPANSION
In the SMOM scheme, the fixed-coupling one-loop dressing functions have explicit analytical expressions given by Eq. (30) where the finite part of the self-energies is written in terms of the adimensional functions F (s), G(s) which are reported in Appendix A.
We can group together the constants in the denominators of Eq. (30) and define the scale-dependent constants so that the dressing functions can be recast as For future reference, we also introduce two differently normalized dressing functions which at s = t satisfy When recast as in Eq.(42), the fixed-coupling dressing functions are known to approach the exact result if the constants F 0 , G 0 take an optimal value [48,49,51]. Only above 2 GeV some deviations occur, as usual in PT, because of the large logarithms ln(p/m) which require a RG improvement of the expansion. It is remarkable that the first derivatives of the inverse dressing functions J 0 (s) −1 = (1 + s −1 ) J(s) −1 and χ(s) −1 , as defined in Eq.(43), do not depend on any parameter and are fully determined, since the constants F 0 , G 0 are canceled by the derivative. Actually, below 2 GeV, all lattice data of different authors, for N = 2, 3, collapse on such curves when the energy units are properly set [49].
By inspection of Eq.(42), we observe that, up to an overall factor, the shape of the dressing functions (as functions of the momentum s) depends on the scaledependent constants F 0 (t), G 0 (t), which is at odds with the concept of multiplicative renormalization. Moreover, since the constants depend on the subtraction point t, they are rather arbitrary and, in principle, even different subtraction points could be taken for the two functions. In previous works, the constants F 0 and G 0 have been used as independent variational parameters for optimizing the one-loop expansion [48,49,51]. The remarkable agreement with the lattice data says that an optimal choice of the constants mimics the neglected higher order terms. Actually, at an optimal subtraction point the higher-order terms, which are also scale dependent, must be so small that the one loop calculation approaches the exact result.
In the SMOM scheme, the constants F 0 , G 0 are not independent since they are defined by Eq.(41) and depend on the scale. It turns out that a very predictive selfcontained variational criterion can be established for determining an optimal subtraction point where the scaledependent constants F 0 (t), G 0 (t) are stationary and the dressing functions approach the exact result. Thus, the optimization of the expansion can be reached by first principles, without any input from the lattice.
If we trust Eq.(42), it is obvious that the ratios J m (t, s)/J m (t ′ , s) and χ m (t, s)/χ m (t ′ , s) in Eq.(32) can be totally independent of the momentum s if and only if the constants F 0 , G 0 are RG invariants i.e.
which is generally not the case for the approximate oneloop expressions of G 0 , F 0 in Eq. (41). Thus, if there is an optimal subtraction point t = t ⋆ where the one-loop dressing functions J m (t ⋆ , s), χ m (t ⋆ , s) approach the exact result, then the functions G 0 , F 0 should be at least stationary at t = t ⋆ , so that the shape of the dressing functions in Eq.(42) does not change for a small shift of the scale around t ⋆ , as required by the multiplicative renormalization of the exact dressing functions. In other words, the optimal scale t ⋆ must be the solution of Eq. (45), which is the point where the constants G 0 , F 0 are locally RG invariant. We expect that any prediction of the one-loop expressions can only be reliable at t = t ⋆ where G 0 (t ⋆ ), F 0 (t ⋆ ) and α(t ⋆ ) get closer to their exact values. An other, related, point of view is that since the shape of the one-loop dressing functions depend on the constants F 0 , G 0 while the shape of the exact dressing functions must be RG invariant, then the best subtraction point must be where the sensitivity to any scale change is minimal, according to Stevenson's principle of minimal sensitivity [65,66]. Again, we find the stationary condition of Eq. (45). The nice thing is that, once we have got a good approximation for the dressing functions and the coupling at t = t ⋆ , we can easily extend the result at any scale t by the exact scaling equations of the SMOM scheme, Eqs. (24), (32), without having to rely on the approximate one-loop relations which are very poor away from t = t ⋆ . From Eq.(32) we can write while from ST identities, Eq.(40) reads yielding an explicit analytical expression for the running coupling in the SMOM scheme. The stationary condition, Eq. (45), can be recast as the solution of two coupled equations for the dressing functions J, χ, as defined in Eq. (43). Using Eq.(44), we observe that the total derivative of J(t, t) and χ(t, t) must coincide, so that if the stationary condition, Eq. (45), is satisfied at t = t ⋆ , then and the scale t ⋆ emerges as the point of tangency of the dressing functions J(t ⋆ , s), χ(t ⋆ , s), while α(t ⋆ ) is the actual value of the functions at the crossing point. The three conditions in Eq. (48) are not enough for determining the four unknown variables t ⋆ , F 0 (t ⋆ ), G 0 (t ⋆ ) and α(t ⋆ ), but they provide a link among them and all of them are fully determined if one is known. We observe that the actual value of the constants F 0 , G 0 depends on the definition of the functions F (s), G(s), since adding a constant to F (or G) and subtracting it from F 0 (or G 0 ), leaves Eqs. (42) and (41) unchanged.
Using the same definitions of Ref. [48] for these functions, as reported in Appendix A, the optimal value F 0 = −0.876 was determined in Ref. [51] by requiring that the complex poles and the whole principal part of the gluon propagator are gauge parameter independent, a feature that the exact gluon propagator inherits from BRST symmetry through the Nielsen identities. For that value, without any fit of parameters, Eqs. (29) and (42) provide a very accurate reproduction of the lattice data of Ref. [18] if the data are scaled in units of m = 0.656 GeV [51]. Since all predictions of the screened expansion are in units of m, taking that value of m is equivalent to share the same energy units of the data sets in Ref. [18].
The ghost dressing function has no poles and no information can be gained by Nielsen identities on the optimal value of G 0 . However, by the stationary condition of Eq.(48), having set F 0 (t ⋆ ) = −0.876 by the analytic properties of the gluon propagator [51], we obtain the optimal set of variables G 0 (t ⋆ ) = 0.14524, t ⋆ = 1.00749 and α(t ⋆ ) = 1.5939 (α s = 4πα/9 = 2.2255 for N = 3).
A graphical solution of the equations is shown in Fig. 1 where the functions J(t ⋆ , s), χ(t ⋆ , s) are plotted together with the running coupling of Eq. (47). As required by Eq.(48), the curves cross at the same point where the dressing functions are also tangent. We observe that the running coupling that arises in the SMOM scheme is monotonic, with a finite IR stable fixed point, at variance with the MOM scheme. It is remarkable that the optimal subtraction point is so close to the mass parameter, µ ⋆ = √ t ⋆ m = 1.004 m. It is not totally unexpected since terms like ln(µ/m) are minimized at that scale.

V. PREDICTIONS OF THE OPTIMIZED EXPANSION
Having found the optimal subtraction point µ 2 /m 2 = t ⋆ for the fixed-coupling screened expansion, we are left with explicit analytical expressions for the propagators and the strong coupling which are fully determined, without any free parameter, from first principles. In this section, the predictions of the optimized expansion are compared with the numerical result of some recent lattice calculations.
In Fig.2 the ghost dressing function χ m (t ⋆ , s) and the gluon propagator m 2 ∆(p) = J m (t ⋆ , s)/(1 + s), renormalized at the optimal scale t ⋆ , are shown together with the large-volume lattice data of Ref. [18] (β = 6.0, L = 80). The lattice data have been multiplicatively renormalized at the scale t ⋆ ≈ 1 and reported in adimensional units by taking m = 0.656 GeV as in Ref. [51]. The plots are not a fit of the data since everything is determined apart from the energy scale m which is obviously dictated by the data. The perfect agreement with the data emerges from the stationary condition, Eq. (45), and from the gaugeparameter independence of the poles, as required by ST identities and Nielsen identities, respectively. All such features are a consequence of the exact BRST symmetry of the total action. Since the symmetry is broken at any finite order of the screened expansion, the one-loop dressing functions give a better approximation at the optimal scale t = t ⋆ where the identities are best satisfied. The prediction of the optimal constant G 0 , from first principles, is encouraging and gives us confidence in the validity of the SMOM scheme for the screened expansion. In fact, if the dressing function were renormalized in the standard MOM scheme, the curves J 0 , χ would not even cross in Fig. 1. The ghost dressing function plays an important role for determining the running coupling by Eq.(47) and the good agreement in Fig. 2 says that the running must be correctly described, at least up to 2 GeV. Beyond that point, the fixed-coupling ghost dressing function is slightly suppressed and some RG improvement is required.
There is an important further prediction that arises from the graphical solution of Eq.(48) in Fig. 1. At the optimal scale t ⋆ the coupling is predicted to be α(t ⋆ ) = 1.5939 which is α s = 2.2255 for N = 3 at µ = µ ⋆ ≈ m. It means that the method is capable of measuring the coupling once the energy units are fixed, as it happens on the lattice. Reversing the argument, if the coupling is known at a given scale, then we get an independent measure of the energy units and m. Thus, it is quite interesting to compare these measures with the lattice data. We observe that, because of Eq.(31) and Eq.(40), we can easily switch from the SMOM to the MOM scheme so that [α s (m)] MOM = 0.5 [α s (m)] SMOM ≈ 1.11. The strong coupling α s , evaluated by Eq.(47) using the predicted value α s (t ⋆ ) = 2.2255 for N = 3, is shown in Fig. 3 together with the lattice data of Ref. [18] for β = 6.0 and L = 80 in the MOM scheme. The data are shown in units of m = 0.656 GeV, as in Fig. 2. In the same figure, using Eq.(49), the calculated strong coupling is also shown in the MOM scheme while the lattice data, which are in the MOM scheme, are also converted and shown in the SMOM scheme. As expected from the agreement in Fig. 2, when rescaled in the MOM scheme, the calculated strong coupling has the same shape of the data, with a maximum occurring at the same point µ/m = 0.85 (µ = 0.56 GeV). However, the estimate of the lattice is 16% higher than the numerical value predicted by the stationary conditions of Eq. (48). As shown in the right panel of Fig. 3, if multiplied by a factor 1.16 the curves collapse on the data. There are good reasons to believe that the disagreement might be caused by the finite lat- tice spacing which is ≈ 0.1 fm for the data set we are considering (β = 6). In fact, as discussed by the same authors of Ref. [18], the numerical value of the coupling depends on the lattice spacing and seems to be overestimated on the lattice, since it decreases with the decrasing of the lattice spacing. Actually, differences larger than 10% have been reported for different lattice spacings. Thus, the prediction of the present calculation is probably closer to the exact result than the lattice esti-mates. At the maximum, the MOM strong coupling is predicted to be α s = 1.157.
For µ > 3m ≈ 2 GeV, the fixed-coupling ghost dressing function is suppressed because of the large logs. In fact, using the explicit expressions of Appendix A, for large values of s, t the function G(x) has the leading behavior G(s) − G(t) ≈ ln(s/t)/4 and the approximate multiplicative renormalization in Eq.(33) is valid for t = t ⋆ ≈ 1 only if α ln s/4 ≪ 1 which is ln s ≪ 2.5 for α(t ⋆ ) ≈ 1.6, yield-ing p ≪ 3.5 m ≈ 2.3 GeV. The problem might be cured by the RG equations, as discussed in Section III. In Fig. 3 we observe a suppression of the running coupling which goes below the data for µ > 3 m, even when rescaled, as in the right panel. The predictions of the fixed-coupling expansion are not reliable beyond that point. We can attempt a matching with the standard running coupling of PT which for a large momentum is very well described by the scheme-independent two-loop beta function and has the approximate solution [70] in powers of 1/ ln µ 2 where, for pure Yang-Mills SU (3) theory, β 0 = 11 and β 1 = 51. As shown in Fig. 3, the approximate two-loop coupling is tangent to the SMOM curve at µ = 2.11 m if Λ = 0.999 m. At that point the strong coupling is α s = 0.59 (g = 2.723). If the curve is rescaled in order to match the data point (right panel), the tangency is found for Λ = 1.08 at µ = 2.31 m where α s = 0.58. Of course, the matching cannot be trusted too much because it occurs at a borderline µ ≈ 2 m which is too large for the fixed-coupling expansion and two small for the two-loop approximation of standard PT. While the data seem to be well described by the approximate two-loop expression of Eq.(51) for µ > 4 m the SMOM curve gives a reliable prediction for µ < 2 m, leaving a narrow window where neither of them can be trusted. That said, it is remarkable that, since Λ ≈ m, the Landau pole is replaced by a mass parameter which gives a scale to the theory and seems to play the same role of Λ.
In the IR, the SMOM scheme seems to be more appealing than the MOM scheme, since the strong coupling saturates instead of decreasing, as one would expect for a phenomenological coupling. As shown in Fig. 3, at variance with the standard MOM scheme, in the screened scheme the strong coupling is a monotonic function of the momentum and reaches its maximum at µ = 0 where α s ≈ 7.03, with a vanishing derivative. Thus, the beta function meets a zero at g = g 0 ≈ 9.40 which is a finite IR-stable fixed point.
Except for the matching range at µ/m ≈ 2−3, Eq.(47) and Eq.(51) provide an explicit analytical expression for the strong coupling and the beta function. No divergence occurs in the IR where no RG improvement seems to be mandatory since the mass parameter screens the theory for µ ≪ m. The beta function is defined as and can be evaluated analytically by a simple derivative.
In the SMOM scheme, the beta is obtained by Eq.(47) and is shown in Fig. 4 for g > 2.723 (µ < 2.11 m) together with the outcome of the approximate two-loop result of Eq.(51) which is shown for g < 2.723. At the matching point, shown as a red circle at g = 2.723, the curves tend to the same value because of the tangency of the couplings in Fig. 3, but a jump in the derivative of the beta is observed. In fact, as discussed before, both curves are not reliable around the matching point and the exact beta function is expected to interpolate smoothly between the two of them. For comparison, in the UV limit g → 0, the exact two-loop beta function of Eq. (50) is also plotted in Fig. 4, showing that the deviation of the approximate Eq.(51) is not small for g > 2.
As shown in Fig. 2 and Fig. 3, the optimized expansion is quite reliable above the matching point, say for g > 3, where the beta reaches a minimum at g ≈ 4.75 and then increases and crosses the zero at the IR stable fixed point g 0 ≈ 9.40. Some numerical values of the beta and of the coupling at specific point locations are listed in Table I.
In Appendix A, explicit expressions are derived for the leading behavior of the functions F (x) and G(x) in the limit x → 0. Deep in the IR, using Eq.(A24), the following explicit expression can be written for the leading behavior of the beta function where g 0 is the IR stable fixed point which can be The leading IR behavior of Eq.(53), valid in the limit g → g0, is shown as a red line. The exact two-loop beta function of Eq. (50), which holds in the UV limit g → 0, is shown as a green line. exactly located as for α(t ⋆ ) = 1.5939 and N = 3, while a, b, c are numerical constants, defined in Appendix A and taking the values The leading behavior of Eq.(53) is also shown in Fig. 4 in the IR limit g → g 0 where it gives a very good approximation of the SMOM beta function for g > 8.5. Asymptotically, the beta tends to the linear behavior β(g) = 2(g − g 0 ).

VI. DISCUSSION
In the IR, the fixed-coupling one-loop screened expansion provides a complete and analytical description of pure Yang-Mills theory if the exact constraints inherited by the BRST symmetry are imposed on the approximate one-loop expressions. While the symmetry is softly broken in the expansion at any finite order, the total action has not been modified, so that the Nielsen and ST identities must be still satisfied by the exact correlation functions of the theory. On the other hand, since the same identities are not satisfied exactly at any finite order, they can be used as a lattice-independent criterion for a choice of the optimal subtraction point.
Since the total action has not been modified, the general renormalization of the theory is not affected by the existence of the mass parameter m which must be regarded as a RG invariant phenomenological scale, fixing the units of the scale-less theory. In that role, the mass parameter m replaces the scale Λ QCD of the standard PT. Actually, the screened expansion and the standard PT share the same general renormalization procedure and the same set of free parameters, since they are different approximations for the same theory. That makes a great difference with other phenomenological models, like the Curci-Ferrari model where a bare mass is added by hand and requires a mass renormalization constant.
The SMOM scheme emerges as the natural renormalization scheme of the screened expansion and the multiplicative renormalization of the theory leads to a stationary condition which gives µ ≈ m as the optimal subtraction point, where the optimal strong coupling is predicted to be α s (m) ≈ 2.23, corresponding to α s (m) ≈ 1.12 in the standard MOM scheme. It is remarkable that, once optimized from first principles, the screened expansion gives an excellent description of the propagators and of the strong coupling if the explicit analytical results are compared with the lattice data in the MOM scheme. Moreover, the lattice-independent measure of the coupling seems to be compatible with the value which can be extrapolated from the lattice data in the continuum limit [18].
The SMOM scheme provides a quite appealing nonperturbative strong coupling which is a monotonic decreasing function of the scale µ/m and saturates in the IR at the finite stable fixed point α s (0) = 7.03 ≈ 2.2 π. An explicit expression for the beta function is provided in the IR with the asymptotic behavior β(g) = 2(g − g 0 ) for g < g 0 . In order to make contact with the phenomenology, it would be interesting to extend the same scheme to the full QCD, including a quark loop in the gluon polarization and comparing the calculated coupling with experimental measures and other approaches [71].
An other important test would arise from a detailed study of the gauge dependence of the coupling in the SMOM scheme. In order to be relevant for the phenomenology, the calculated strong coupling should be gauge invariant, but that would be hard to tell by a oneloop calculation. The Landau gauge is quite special since Z c 1 = 1, exactly, because of the non-renormalization theorem. In a generic covariant gauge, the SMOM scheme would require an explicit calculation of the vertex for determining the finite part of Z c 1 , which would also depend on the gauge parameter ξ. While the total gauge dependence might not cancel exactly in the coupling, the residual dependence might not be genuine as it could arise from the further one-loop approximation for the vertex. In that respect, the Landau gauge provides our best estimate for the strong coupling. where the logarithmic functions L x are and the rational parts are where the logarithmic functions L ′ x , for x = a, b, c, are L ′ a (x) = 6x 4 − 16x 3 − 68x 2 + 80x + 144 and R ′ (x) is the sum of all the rational terms coming out from the derivatives R ′ (x) = 12 x + 106 Observe that the functions L ′ x , R ′ are not the derivatives of the corresponding functions L x , R.

IR leading behavior
Deep in the IR, the following leading behavior is obtained in the limit x → 0 which has been checked to be a very good approximation for x < 0.01 (p < 0.1 m). At the same order of approximation, the dressing function J(t ⋆ , x) can be written as which holds for x < 0.015, with the constants: The leading behavior of the function L gh (x) is yielding the approximate expression which has been checked to be a very good approximation for x < 0.2, with the constants k 0 = 5 24 + G 0 , k 1 = 2 9 , k L = 1 6 . (A13) At the same order of approximation, the dressing function χ(t ⋆ , x) can be written as which holds for x < 0.01, with the constants: By insertion in Eq. (47), deep in the IR, the leading behavior of the SMOM running coupling is and g 0 is the IR stable fixed point for α(t ⋆ ) = 1.5939 and N = 3.
The leading behavior of the beta function in the IR follows by its definition which by Eq.(A16) gives the coupled equations β(g) = 2(g − g 0 ) − 2b g 0 x(g) x(g) = f −1 g − g 0 g 0 (A20) where f (x) = x (a− b ln x). At the same order of approximation, the function f (x) can be inverted by iteration. Denoting by ǫ the small positive variable and observing that f [x(g)] = b ǫ, the function f (x) can be written as where η = a/b. Then, iterating, x(g) = ǫ η − ln ǫ + ln(− ln x) + ln 1 − η ln x