Analytical study of Yang-Mills theory in the infrared from first principles

Pure Yang-Mills SU(N) theory is studied in the Landau gauge and four dimensional space. While leaving the original Lagrangian unmodified, a double perturbative expansion is devised, based on a massive free-particle propagator. In dimensional regularization, all diverging mass terms cancel exactly in the double expansion, without the need to include mass counterterms that would spoil the symmetry of the Lagrangian. No free parameters are included that were not in the original theory, yielding a fully analytical approach from first principles. The expansion is safe in the infrared and is equivalent to the standard perturbation theory in the UV. At one-loop, explicit analytical expressions are given for the propagators and the running coupling and are found in excellent agreement with the data of lattice simulations. A universal scaling property is predicted for the inverse propagators and shown to be satisfied by the lattice data. Higher loops are found to be negligible in the infrared below 300 MeV where the coupling becomes small and the one-loop approximation is under full control.


I. INTRODUCTION
In modern textbooks on QCD, the infrared domain is usually called non-perturbative just because standard perturbation theory breaks down at the low-energy scale Λ QCD ≈ 200 MeV. While the high energy behaviour of the theory is under control and an analytical study of non-Abelian gauge theories is usually achieved by a loop expansion in the UV, no analytical first-principle description of the infrared can be found in books where the subject is usually discussed by phenomenological models that rely on numerical lattice simulations.
It is now widely believed that in the Landau gauge the gluon propagator is finite and an effective coupling can be defined that is infrared safe and relatively small. As discussed by Cornwall [30] in 1982, the gluon may acquire a dynamical mass in the infrared without breaking the gauge invariance of the theory. The effect cannot be described by the standard perturbation theory at any finite order because of gauge invariance that makes the polarization transverse and prohibits any shift of the pole in the gluon propagator. That is one of the reasons why the standard perturbation theory cannot predict the correct phenomenology in the infrared.
An other reason is the occurrence of a Landau pole in the running of the coupling that makes evident the failure of the perturbative expansion below Λ QCD . However, in the Landau gauge the ghost-gluon vertex function can be shown to be finite [29] and a running coupling can be defined by the product of two-point correlators. Being massive, the gluon propagator is finite and its dressing function vanishes in the infrared yielding a finite running coupling that reaches a maximum and decreases in the low-energy limit [24]. On the other hand, if the Landau pole is an artifact of the perturbative expansion, the relatively small value of the real coupling suggests that we could manage to set up a different perturbative scheme in the infrared. Actually, in order to make physical sense, a perturbative expansion requires that the lowest order term should approximately describe the exact result. While that condition is fulfilled by the standard perturbation theory in the UV, where the propagator is not massive, the dynamical mass of the gluon makes the free propagator unsuitable for describing the low energy limit. Thus we would expect that, by a change of the expansion point, a perturbative approach to QCD in the infrared could be viable.
There is some evidence that inclusion of a mass by hand in the Lagrangian gives a phenomenolgical model that describes very well the lattice data in the infrared at one loop [31][32][33]. However that model could be hardly justified by first principles because of the mass that breaks the gauge invariance of the Lagrangian. Even in a fixed gauge, BRST symmetry is broken and a mass counterterm must be included for renormalizing the theory, thus introducing spurious free parameters in the model.
A change of the expansion point can be achieved by first principles without changing the original Lagrangian. Variational calculations have been proposed [16][17][18][19] where the zeroth order propagator is a trial unknown function to be determined by some set of stationary conditions. The added propagator is subtracted in the interaction, leaving the total action unchanged. The idea is not new and goes back to the works on the Gaussian effective potential [34][35][36][37][38][39][40][41][42][43][44][45][46] where an unknown mass parameter was inserted in the zeroth order propagator and subtracted from the interaction, yielding a pure variational approximation with the mass that acts as a variational parameter. Some recent variational calculations on Yang-Mills theory [18,19] have shown that, provided that we change the expansion point, a fair agreement with the lattice data can be achieved without too much numerical effort. Thus we expect that it is not very important the actual choice of the zeroth order propagator provided that it is massive. A simple free-particle Yukawa-type massive propagator would be enough and the corrections due to the interaction would then be manageable by perturbation theory.
A first attempt along these lines was reported in Ref. [20] where, by a second order massive expansion, the gluon and ghost propagators are evaluated and found in fair agreement with the lattice data. The integrals were regularized by a simple cutoff that breaks the BRST symmetry and gives rise to several drawbacks like quadratic divergences and the need of a mass counterterm. However, a fine tuning of the mass parameter seems to cure the drawbacks yielding an optimized expansion that reproduces the lattice data.
In this paper, the difficulties of dealing with a cutoff are avoided by the use of a more robust dimensional regularization scheme, yielding a more rigorous perturbative study of pure SU (N ) Yang-Mills theory from first principles. While the original Lagrangian is not changed in any way, the outcome is a one-loop analytical description that is infrared safe and in striking agreement with the data of lattice simulations. Moreover the result can be improved by including higher-order terms and by use of standard Renormalization Group (RG) techniques for reducing the effect of higher order terms.
A very interesting property of the massive expansion is the cancellation of all diverging mass terms without including any spurious mass counterterm. Only wave function renormalization constants are required and, in the minimal subtraction scheme, these constants are the same of the standard perturbative expansion, thus ensuring that the correct UV behaviour is recovered.
The massive expansion is discussed in the Landau gauge in the present paper. The Landau gauge is probably the optimal choice for the expansion, because of the transversality of the propagator that makes the longitudinal polarization irrelevant. In the Landau gauge the problem decouples and a fully analytical result can be found for the propagators at one-loop.
While massive models have been studied before and found in good agreement with the data of lattice simulations [31][32][33], the present calculation is very different because the Lagrangian is not modified, overall BRST symmetry is not broken and no free parameters are added to the exact Yang-Mills theory, yielding a description that is based on first principles and can be improved order by order. Thus, at variance with previous massive models, the present method would not give a mass to the photon. The paper is organized as follows: in Section II the massive expansion is developed for pure SU (N ) Yang-Mills theory in a generic covariant gauge; in Section III the double expansion is set up in the Landau gauge; in Section IV the explicit cancellation of the diverging mass terms is discussed in detail; in Section V explicit analytical expressions are derived for the propagators at oneloop; in Section VI the one-loop propagators and their scaling properties are compared with the available lattice data; in Section VII the running coupling is evaluated and its sensitivity to the renormalization conditions is discussed, showing that the approximation is under full control below 300 MeV; finally, in Section VIII the main results are discussed. Explicit analytical expressions for the propagators are given in the Appendix.

II. SET UP OF THE MASSIVE EXPANSION IN A GENERIC GAUGE
Let us consider pure Yang-Mills SU (N ) gauge theory without external fermions in a d-dimensional space. The Lagrangian can be written as where L Y M is the Yang-Mills term L f ix is a gauge fixing term and L F P is the ghost term arising from the Faddev-Popov determinant.
In terms of the gauge fields, the tensor operatorF µν iŝ and the generators of SU (N ) satisfy the algebra with the structure constants normalized according to If a generic covariant gauge-fixing term is chosen the total action can be written as S tot = S 0 + S I where the free-particle term is and the interaction is with the three local interaction terms that read In Eq. (8), ∆ 0 and G 0 are the standard free-particle propagators for gluons and ghosts and their Fourier transforms are Here the transverse and longitudinal projectors are defined as where η µν is the metric tensor. A shift of the pole in the gluon propagator can be introduced by an unconventional splitting of the total action. Since we have the freedom of adding and subtracting the same arbitrary term δS to the total action we can take where the vertex function δΓ is a shift of the inverse propagator and ∆ m µν is a massive free-particle propagator Here the dynamical mass M(p) and the longitudinal function A(p) are left totally arbitrary. While the total action cannot depend on them, just because δS is added and subtracted again, any expansion in powers of the new shifted interaction S I → S I − δS is going to depend on the choice of δS because of the approximation. Thus, it is the approximation that changes but we are not changing the content of the exact theory. The shift δS has two effects: the free-particle propagator ∆ 0 µν is replaced by the massive propagator ∆ m µν in S 0 ; a counterterm −δS is added to the interaction S I .
From now on, let us drop all color indices in the diagonal matrices. Inserting Eq. (11) and (16) in Eq.(15) the counterterm reads δΓ µν (p) = M(p) 2 t µν (p) + A(p)ℓ µν (p) (17) and must be added to the standard vertices arising from Eq.(10). The proper gluon polarization Π and ghost self energy Σ can be evaluated, order by order, summing up Feynman graphs where ∆ m µν is the zeroth order gluon propagator. By Lorentz invariance we can write and the dressed propagators are where the transverse and longitudinal parts read At tree level, the polarization is just given by the counterterm δΓ of Eq. (17), so that the tree-terms Π T tree = M 2 , Π L tree = A just cancel the shifts in the dressed propagator ∆ of Eq. (20), giving back the standard free-particle propagator of Eq. (11). In fact, at tree-level, nothing really changes.
Summing up all loops, the exact dressed propagator can be written as As a consequence of gauge invariance, the exact longitudinal polarization Π L loops must be zero and the longitudinal part of the exact propagator must be equal to its tree-level value ∆ L = −ξ/p 2 , just because the loopterms cannot change it, as recently confirmed by lattice simulations [27]. Since Π L loops and Π T loops are evaluated by insertion of the modified propagator ∆ m µν in the loops, they can be considered as functionals of the arbitrary functions M, A. Thus, summing up all loops, the following constraints must hold for the exact polarization functions: Expanding in powers of the total interaction, including the counterterm δΓ among the vertices and writing down the Feynman graphs, we can truncate the expansion at a finite order yielding approximate functionals that may not satisfy the constraints of Eq. (22) exactly. For instance, the exact vanishing of the transverse polarization would be lost unless BRST symmetry is mantained order by order. Actually, while the total Lagrangian has not been changed and mantains its symmetry, the two parts S 0 and S I might be not BRST invariant because of the arbitrary shift δS. Then the exact symmetry is lost in the expansion at any finite order and the constraints are expected to hold only approximately unless all the graphs are summed up. The outcome of the truncated expansion becomes sensitive to the choice of the functions A, M thus suggesting the use of Eq.(22) as variational stationary conditions.
Since we expect that the approximation should work better if the zeroth order propagator ∆ µν m is a good approximation of the exact one, then comparing Eq. (21) and Eq.(16), a self-consistent method could be set up by requiring that Summing up all the loops these equations would be equivalent to SDE. Variational methods of this kind have been investigated in several works [16][17][18][19] and require the solution of integral equations that can hardly be treated analytically.
In the Landau gauge the problem decouples and a fully analytical result can be found for the propagators. The Landau gauge is very special as the gluon propagator is transverse and does not depend on the choice of the function A. In the limit ξ → 0 the longitudinal part ∆ L is exactly zero in Eq. (21) and is decoupled from the longitudinal polarization that becomes irrelevant for the calculation of the propagator. In fact, in the same limit, the zeroth order propagator ∆ m µν in Eq.(16) becomes transverse for any choice of A and the longitudinal part of the counterterm δΓ does not give any contribution in the loops when sandwiched by two transverse propagators. Thus the only action of A is at tree level where it cancels itself in the gluon propagator and disappears. Then, in the Landau gauge, we can safely drop all longitudinal parts and set A = 0 without affecting the calculation. Because of the decoupling, the calculation of the longitudinal and of the transverse parts can be seen as two separate problems that may even require different orders of approximation. That simplifies things considerably, since a poor approximation for the longitudinal polarization would not affect the accuracy of the propagator. That also explains why reliable results for the propagator can be achieved even when the BRST symmetry is broken [32] in the Landau gauge. Moreover, since the total Lagrangian has not been modified, the overall BRST symmetry is unbroken and the constraints in Eq.(22) must be satisfied asymptotically. Thus, if required, a better approximation for the longitudinal polarization can always be achieved in a separate calculation by adding more terms to the expansion.

III. DOUBLE EXPANSION IN THE LANDAU GAUGE
The Landau gauge is probably the optimal choice for the massive expansion, as discussed in the previous section. In the limit ξ → 0 the gluon propagators are transverse exactly and, having set A = 0, we can simplify the notation and drop the projectors t µν everywhere whenever each term is transverse. Moreover, we make the minimal assumption of taking the arbitrary function M equal to a constant mass scale M = m. In fact, variational calculations seem to suggest [18,19] that the actual form of the zeroth order propagator is not important provided that it is massive. A constant mass simplifies the calculation and allows the use of dimensional regularization.
We can use the standard formalism of Feynman graphs with a massive zeroth order propagator that reads and a counterterm that must be added to the standard three-particle ghostgluon and gluon-gluon vertices of order O(g) and to the four-particle gluon-gluon vertex of order O(g 2 ) according to Eq.(10). Thus, the total interaction is a mixture of terms that depend on the coupling strength g and a counterterm that does not vanish in the limit g → 0.
A perturbative expansion in powers of the total interaction would contain at any order different powers of g but the same number of vertices (including the counterterm among vertices) and we may define the order of a term as the number of vertices in the graph. Of course, we could easily sum up some infinite set of graphs, like the chain graphs in Fig.1. Formally, summing up all graphs with n insertions of the counterterm in the internal gluon lines, would cancel the pole shift and would give back the standard perturbation theory In fact, this is what we would get exactly at tree level, as discussed in the previous section. That just says that the massive and the standard expansions are equivalent if we sum up all graphs. On the other hand, at any finite order, the massive expansion is not equivalent to  the standard perturbation theory, but the two expansions differ by an infinite class of graphs that amounts to some non-perturbative content. For instance, the massive zeroth order propagator ∆ m cannot be obtained by the standard perturbation theory at any finite order because of the gauge invariance of the theory that does not allow any shift of the pole. We may reverse the argument and observe that, while in the UV the geometric expansion in Eq.(26) is convergent and the two perturbation theories must give the same result, when p 2 → m 2 each single term in Eq.(26) diverges and the formal sum of infinite poles amounts to some non-perturbative content that makes the theories different at any finite order. We can predict that the scale m should be close to the Landau pole Λ where the standard perturbation theory breaks down.
Since we know that the gluon develops a dynamical mass in the infrared, we do not want to sum the chain graphs in Fig.1 but prefer to truncate the power expansion at some finite order. An expansion in powers of the total interaction S I is more efficient than the standard expansion in powers of the coupling g. The counterterm δΓ has the important effect of reducing the weight of the total interaction since, in principle, if the zeroth order propagator were exact, the total polarization would be exactly zero. For that reason, we define the order of a graph as the number of vertices that are included, reflecting the power of S I rather than the number of loops. Thus the tree-level graphs must be regarded as first order. As shown in Fig.2 the one-loop tadpole (1b) is first order while the gluon loop (2b) is second order. Any insertion of the counterterm δΓ increases the order by one.
If the effective coupling is small, as it turns out to be according to non-perturbative calculations, not all the graphs have the same weight in the expansion. Since the number of loops is equal to the power of g 2 in a graph, two-loop graphs must be much smaller than one-loop and tree graphs. We can consider a double expansion in powers of the total interaction and in powers of the coupling: we expand up to the nth order, retaining graphs with n vertices at most, and then neglect all graphs with more Figure 2: Two-point graphs with no more than three vertices and no more than one loop. In the next sections, the ghost self energy and the gluon polarization are obtained by the sum of all the graphs in the figure. than ℓ loops. In Fig.2 the lower order graphs are shown up to n = 3 and ℓ = 1.
A very important feature of the double expansion is that there is no need to include mass counterterms for regularizing the divergences. All diverging mass terms cancel exactly in the expansion. Thus we avoid to insert spurious mass parameters that were not in the original Lagrangian. The cancellation can be easily explained by the following argument. Since we did not change the Lagrangian and it was renormalizable (without mass counterterms because of BRST symmetry), all the diverging mass terms must cancel if we sum up all graphs. In fact, no diverging mass term arises in the standard perturbative expansion. The cancellation must be given by the sum of infinite graphs with counterterm insertions δΓ in the loops that, according to Eq.(26) and Fig.1, restore the pole of the propagator and cancel the mass. However, if we inspect the graphs in Fig.2, we can easily see that any insertion of δΓ in a loop reduces the degree of divergence of the graph so that they become finite after a finite number of insertions. Thus, if the divergences must cancel, they will cancel at a finite order of the expansion provided that we retain more counterterm insertions than loops. If n is large enough, then all divergences in the mass terms are cancelled by the counterterms in the loops. For instance, at one loop we only need n = 3 as shown in Fig.2.
While in this paper we report the results for a oneloop (third-order) approximation, the extension to higher loops is straightforward and the regularization follows the same path of standard perturbation theory, with all the divergences that can be cancelled by the usual wave function renormalization constants.

IV. CANCELLATION OF MASS DIVERGENCES IN DIMENSIONAL REGULARIZATION
The exact cancellation of the diverging mass-terms can be carried out explicitely in dimensional regularization expanding in powers of ǫ = 4 − d.
The insertion of one counterterm in a loop can be seen as the replacement in the internal gluon line. If there are no other counterterm insertions in the same graph, then the dependence on m 2 must come from the massive propagators and a derivative of the whole nth-order ℓ-loop graph gives the sum of all (n+ 1)-order ℓ-loop graphs that can be written by a single insertion of δΓ in any position.
In dimensional regularization, any diverging mass term that arises from a loop can be expressed as a pole c m 2 /ǫ where c is a factor. Inserting this term in Eq. (27), we see that a counterterm in the loop gives a crossed-loop graph with the opposite diverging term −cm 2 /ǫ. The argument also suggests a simple way to evaluate the crossed-loop graphs by Eq. (27).
At one-loop, we must truncate the expansion at the order n = 3 for a full cancellation of all the diverging mass terms. While higher-order terms could be included without introducing any further divergence at one-loop, in this paper we explore the minimal approximation and sum up all graphs up to n = 3 as shown in Fig.2. It is not difficult to show that in the limit p → 0 the gluon polarization is finite but not zero. The existence of a finite limit Π(0) = 0 is crucial for the existence of a finite gluon propagator in the infrared.
First of all, let us evaluate the constant graphs at the order n = 3. At the lowest order (n = 1, ℓ = 0) the counterterm δΓ gives the constant graph Π 1a = m 2 that cancels the shift of the pole in the propagator. Exact integral expressions for the loop graphs have been reported by other authors in the Landau gauge. In Ref. [19] all one-loop graphs are reported for any gauge, any space dimension and any choice of the zeroth order propagator. In Landau gauge and Euclidean space, the constant tadpole Π 1b can be written as Expanding around d = 4, in the M S scheme, having hided the factor N inside an effective coupling α defined as The crossed tadpole Π 1c follows by a derivative according to Eq.(27) As expected, the diverging terms cancel in the sum Π 1b + Π 1c . In fact, the double-crossed tadpole Π 1d is finite and including its symmetry factor it reads so that the sum of the constant graphs is While the ghost loop vanishes in the limit p → 0, a finite Π(0) = 0 can also arise from the gluon loop Π 2b that in the Landau gauge (in Euclidean space) can be written as [19] where k 2 ⊥ = [k 2 − (k · p) 2 /p 2 ] and the kernel F can be decomposed as The calculation of this graph is straightforward but tedious. The integral can be evaluated analytically and the result is reported in the next section. If we take the limit p → 0 before integrating, we find that F → 1 and a mass term arises The integral is trivial and expanding around d = 4 in the M S scheme we can write it as Adding the crossed loop Π 2c with its symmetry factor, the divergences cancel (38) and adding the constant graphs in Eq.(33), the one-loop dressed propagators can be written as While the explicit calculation requires the evaluation of the gluon and ghost loop and of the ghost self-energy, we observe that a finite mass-term has survived in the cancellation, so that the dressed propagator ∆(0) −1 = 5αm 2 /8 is finite and of order α. Actually, we checked that the full propagators in Eq.(39), when renormalized, do not depend on the precise value of the factor 5/8 that arises by truncating the expansion at the third order. A minor change of that coefficient is absorbed by a change of the mass parameter and of the renormalization constants without affecting the final result. That is an important feature since otherwise the whole calculation would depend on the somehow arbitrary truncation of the expansion. In fact, while higher-order terms would add very small corrections in the UV because of the factor (−p 2 + m 2 ) n+1 ∼ (−p 2 ) n+1 in the denominators of Eq.(26), in the limit p → 0 the corrections might not be negligible. In that limit we find a hierarchy in the significance of the crossed terms. The most important effect arises at tree-level since the treegraph Π (1a) in Fig.2 cancels the entire shift of the pole in the propagator, as discussed in Section II. Thus, a finite Π(0) = 0 can only arise from loops and the massive expansion would not predict any mass for the photon. At one-loop, a first insertion of the counterterm gives diverging crossed graphs that cancel the divergence of the loops entirely. Inclusion of those terms is crucial for the renormalization of the theory. On the other hand, the insertion of n counterterms in a loop, with n ≥ 2, gives finite terms that only add some fractions of αm 2 to Π(0). These terms decrease as ∼ 1/n 2 and subtract each other, with a positive series of terms coming from the tadpole graph and a negative series arising from the gluon loop. Thus, the inclusion of higher order terms would only give a slight decrease of the coefficient of Π(0) = −5αm 2 /8 in Eq. (39). That change is compensated by an increase of the mass parameter and by a change of the renormalization constants, without making any real difference in the renormalized propagators. In that sense, the minimal choice of a third order expansion has nothing special in itself and no dramatic effect is expected if higher-order terms are included.

V. ONE-LOOP PROPAGATORS
The explicit evaluation of the propagators at one-loop and order n = 3 requires the sum of the gluon loop Π 2b , the crossed loop Π 2c and the ghost loop Π 2a for the gluon propagator and the sum of the one-loop and crossed-loop self energy graphs for the ghost propagator, as shown in Fig.2.
From now on, we switch to Euclidean space, expand the graphs around d = 4 in the M S scheme and use the adimensional variable s = p 2 /m 2 . The gluon loop Π 2b is given by the integral in Eq.(34) that gives a diverging part and a finite part where L A , L B are the logarithmic functions L A (s) = (s 2 − 20s + 12) 4 + s s The ghost loop Π 2a is a standard graph and in the Landau gauge it is given by the integral [19] The integral is straightforward and the diverging part is while the finite part reads The first of self-energy graphs in Fig.2, the standard one-loop graph, is given by the integral [19] that yields a diverging term and a finite part where the function g(s) is If we do not add the crossed loops and take the sum of the finite parts, Eqs. (41) and (45), then we recover the finite part of the one-loop polarization function where the function f (s) is We observe that in the limit s → 0 the logarithmic functions have the limits sL A (s) → −96 and s(L B (s) − 2s −2 ) → −15, so that sf (s) → −111.
The crossed loops can be included very easily by a derivative with respect to m 2 , as discussed in the previous section below Eq.(27) where we include all finite and diverging parts in the derivative. The derivative of the diverging parts gives the finite terms −αp 2 /4 and 13αp 2 /18 that must be added to the finite parts of self-energy and polarization function, respectively. Thus the diverging parts do not change and are given by the one-loop terms, Eqs. (47) and (52).
Performing the derivative of the finite parts we obtain where f ′ and g ′ are the derivatives of f and g, respectively. The bare propagators follow by the insertion of finite and diverging parts in Eq. (39). The propagators can be made finite by the standard wave function renormalization. At one loop, the only residual mass term is finite and of order α, so that the divergences in Eq.(39) are absorbed by the wave function renormalization constants Z A , Z ω . In the M S scheme we find by Eqs. (47) and (52) thus reproducing the same UV behaviour of the standard one-loop approximation.
It is useful to introduce the adimensional ghost and gluon dressing functions that, once renormalized by the constants Z A , Z ω , are finite and read where F (s) = 5 8s while in the infrared, for s → 0, we find that G(s) tends to a constant and F (s) ≈ 5/(8s), so that χ(0) is finite and J(s) ≈ 8s/(5α), yielding ∆(0) −1 = 5αm 2 /8 as expected from Eq. (39). We observe that in the UV, the asymptotic behaviour of Eq.(59) is precisely what we need for canceling the dependence on m in the dressing functions. In fact, in the UV, Eq.(57) can be written as which is the standard UV behaviour that we expected by inspection of the renormalization constants Eq.(55). The constants in Eq.(57) have no direct physical meaning and depend on the special choice of renormalization constants in the M S scheme. We can subtract the dressing functions at a generic point s 0 and, without fixing any special renormalization condition, we can write them in the more general form that extends the standard UV one-loop behaviour of the Eqs.(60), sharing with them the same asymptotic behaviour for s, s 0 ≫ 1 according to Eq.(59). We observe that in general, we might not have the freedom of setting J(s 0 ) = χ(s 0 ) = 1 in Eq.(61). Actually, F (s) is not a monotonic function, it has a minimum and is bounded from below, so that J(s) −1 must also be bounded in Eq.(61). Of course, that is just a limit of the one-loop approximation and the dressing functions can be renormalized at will by a different choice of the renormalization constants. The point is that if the dressing functions are multiplied by the arbitrary factor Z = 1 + αδZ then, at one-loop, that is equivalent to the subtraction of αδZ on the right-hand sides of the Eqs.(57). That only makes sense if δZ is small and Z ≈ 1. While in principle Z can take any value, even much larger or smaller than 1, the one-loop subtraction can only compensate a small value of δZ. That is not a problem in Eq.(61) provided that we take account of any large renormalization factor by direct multiplicative renormalization of χ(s 0 ) and J(s 0 ). Then, if the energy s is not too far from the subtraction point s 0 , the one-loop correction is small as it must be. An important consequence is that by Eq.(61) we can predict that at one-loop, up to an arbitrary multiplicative renormalization constant, the inverse dressing functions are given by the universal functions F (s) and G(s) up to an additive renormalization constant. Such scaling property is satisfied quite well by the lattice data, thus enforcing the idea that perturbation theory can provide important insights on QCD in the infrared.

VI. SCALING PROPERTIES ON THE LATTICE
The predictive content of the theory can be tested by a direct comparison with the lattice data. First of all, we would like to explore the scaling properties that emerge from Eq.(61) and that seem to be satisfied by the available lattice data for SU (2) and SU (3). In fact, in Eq.(61) any dependence on α is absorbed by the multiplicative renormalization constants of χ and J. By such renormalization, the inverse dressing functions are entirely determined by the universal functions F (s) and G(s) up to an additive constant. In other words, by a special choice of the renormalization constants, all dressing functions can be translated on top of the same curve by a vertical shift. In order to make that more explicit, we can write Eq.(61) as where F 0 and G 0 are a pair of constants depending on the subtraction point s 0 , on the bare coupling and on the normalization of the dressing functions χ(s 0 ), J(s 0 ), while Z G , Z F are arbitrary renormalization constants that also absorb the dependence on α. While these equations predict a scaling property that is a stringent test for the oneloop approximation, the predictive content is remarkable: the derivatives of the inverse dressing functions must be equal to the derivatives of the universal functions F (s), G(s) up to an irrelevant multiplicative factor, while the additive constants F 0 , G 0 emerge as unknown integration constants. The mass parameter m provides the natural energy units that cannot be predicted by the theory and can only be fixed by comparison with physical observables or lattice data. In fact, the total Lagrangian does not contain any energy scale and, as for lattice calculations, the natural scale must be regarded as a phenomenological quantity. However, once the mass m is fixed, the original arbitrariness of its choice is reflected in a spurious dependence on the subtraction point s 0 which is the only scale that remains free in the theory. We expect that the residual dependence on s 0 , which is implicit in the constants F 0 , G 0 , should decrease if the approximation  Table I  is improved by the inclusion of higher loops.  The function F (s) + F 0 is shown in Fig.3 together with the lattice data for the gluon inverse dressing function. For SU (3) the data points are extracted from a figure of Ref. [24] while for SU (2) the interpolation function of Ref. [47] is used, valid in the range 0.7-3.0 GeV. The data are scaled by the renormalization constants in Table I and shown to collapse on the one-loop function F (s) by a vertical translation. Eq.(62) is satisfied very well in the whole range of the lattice data. There is a pronounced minimum that fixes the energy scale at m = 0.73 GeV for SU (3) and m = 0.77 GeV for SU (2). In Fig.3 and Fig.4, the energy units of the data for SU (2) have been scaled by the ratio of the masses in order to superimpose them on the data for SU (3). An enlargement of the area of the minimum is shown in Fig.4 where the deviations between the curves are amplified but found to be smaller than the fluctuations of the lattice data.
The function G(s)+G 0 is shown in Fig.5 together with the lattice data for the ghost inverse dressing function. As in Fig.3, the lattice data for SU (3) are extracted from a figure of Ref. [24] while the data for SU (2) are given by the interpolation function of Ref. [48], valid in the range 0.2-3.5 GeV. Again, the data are scaled by the renormalization constants in Table I and collapse on the one-loop function G(s) by a vertical translation. The energy units are the same of Fig.3 and Fig.4, i.e. the same values of m are required for ghost and gluon dressing functions. We can see that the scaling properties predicted by Eq.(62) are also satisfied very well by the ghost dressing function.
Overall, we find a very satisfactory description of the lattice data if the renormalization constants, the multiplicative factors Z F , Z G and the additive constants F 0 , G 0 are fixed as in Table I. While the multiplicative fac-tors are not relevant anyway, we find a slight dependence on the additive constants that cannot be compensated by a change of the factors, because of the one-loop approximation. Once the energy scale m is fixed, no other free parameters are left besides the renormalization constants, so that the agreement with the lattice data is remarkable and really encouraging.
The gluon propagator and the ghost dressing function seem to show an even better accuracy than their inverse, because of the scale. For instance, for SU (3) the gluon propagator and the ghost dressing function are reported in Fig.6 and Fig.7, respectively, together with the lattice data of Ref. [24]. The renormalization constants are set at the same values of Table I as discussed above. We observe that the gluon propagator is not convex. Actually, it is not even a monotonic function of p, as shown in Fig.8 where an enlargement of the deep infrared area is displayed in more detail. That property is usually assumed to be a sign of confinement. A comparison of the gluon propagator with the lattice data of Ref. [47] for SU (2) is given in Fig.9.

VII. RUNNING COUPLING
In the Landau gauge the ghost-gluon vertex is regular [29] and the vertex renormalization constant can be set to one in a momentum-subtraction scheme, so that a runnig coupling is usually defined by the RG invariant product of the dressing functions Having reproduced the dressing functions very well, we expect a very good agreement with the lattice for the (64) and depends on the renormalization point µ = µ 0 where we set α s (µ 0 ) at a given phenomenological value. We can renormalize the coupling at the point µ = 2 GeV where the lattice data of Ref. [24] give α s = 0.37 for SU (3). That is a good compromise as the coupling is still quite small while the energy is not too large, so that we can still neglect the RG effects that become important in the UV limit [32]. We will refer to this point as the large energy renormalization point. Using the values of Table  I for F 0 , G 0 and m at N = 3, the running coupling of Eq.(64) is displayed in Fig. 10, together with the lattice data of Ref. [24]. The agreement is very good in the whole infrared range for µ < 2.5 GeV. In the UV, when µ > 2.5 GeV, we observe that Eq.(62) starts to deviate from the lattice data. That is a known problem that can be cured by a consistent running of the coupling in the one-loop calculation according to the RG equations, as shown in Ref. [32]. On the other hand, in the infrared the agreement is impressive for a one-loop calculation.
It is instructive to explore how sensitive the result is to the choice of the additive renormalization constants F 0 , G 0 , which are the only free parameters of the calculation. From a physical point of view, we would expect that if the running coupling α s (µ) is the true effective coupling at the scale µ, then the one-loop approximation should be working very well deep in the infrared where α s → 0. That would be very interesting for future perturbative work. A test of the one-loop approximation comes from the sensitivity to changes of the additive constants. If the approximation is under full control, then any small change of F 0 and G 0 should be compensated by the multiplicative renormalization constants, thus canceling in the normalized ratio of Eq.(64). In Fig.11, the grey pattern shows the area spanned by the running coupling α s (µ) of Eq.(64) when the additive renormalization constants F 0 , G 0 are changed by ±25% around the values of Table I. Ignoring RG effects in the UV and comparing with the best running coupling of Fig.10, which is also  Figure 11: The filled grey pattern is the area spanned by the coupling αs(µ) when the constants F0, G0 are changed by ±25% with respect to the values in Table I  shown in the figure, we see that the deviations are very small in the UV and start growing up when α s ≈ 0.6. They increase until α s reaches its maximum and then decrease getting smaller and smaller in the infrared limit µ → 0. That enforces the idea that, deep in the infrared, the one-loop approximation could be under full control. Moreover, the sensitivity to the additive constants seems to be even smaller in the infrared if the renormalization point is taken at a very low energy. In Fig. 12 the deviations are evaluated as before, by Eq.(64), but renormalizing the coupling at µ = 0.15 GeV where α s = 0.2. We can see that the running coupling seems to be not sensitive at all to the choice of the additive constants until α s ≈ 0.6, and the approximation seems to be under full control below 300 MeV. In other words, regardless of the actual value of the renormalization constants, all curves evaluated by Eq.(64) collapse on the lattice data below 300 MeV. That is a remarkable feature as, by a proper choice of the renormalization point, the present one-loop approximation provides a very accurate description of the running coupling below 300 MeV (Fig. 12) or above 1.5 GeV (Fig. 11), without adjusting any free parameter, from first principles. In the range between 0.3 and 1.3 GeV, where α s > 0.6, Eq.(64) can be still tuned on the lattice data, as shown in Fig. 10, but the increased sensitivity to the additive renormalization constants is a sign of the limits of the one-loop approximation. However, since the calculation is from first principles, we expect that the sensitivity to the additive constants should decrease when higher loops are included in the expansion. From a technical point of view, Eq.(64) provides a very good interpolation of the lattice data and is not sensitive to the choice of the renormalization constants below 300 MeV and above 1.5 GeV. Then, it could make sense to introduce a third fixed point by just renormalizing at the scale where the deviations are larger and pinpoint α s at its maximum. If we renormalize at the maximum point µ = 0.67 GeV setting α s = 1.21, the deviations are quite small over the whole range of energies, as shown in Fig.13. That suggests that we can get rid somehow of the additive constants and write some universal function for the running coupling, free of any parameter, albeit slightly approximate.
Let us pretend that we can set χ(s 0 ) = J(s 0 ) = 1 in Eq.(61) and insert it in Eq.(63). Then, neglecting higher powers of α, the running coupling takes the simple shape where α(p 2 /m 2 ) = 3N α s (p)/(4π) and the universal function S(s) is defined as and does not contain any free parameter. In the UV, the running coupling α(s) incorporates the standard one-loop leading behaviour. In fact, by Eq.(59), for s, s 0 ≫ 1 which does not depend on the scale m. In the infrared, the function S(s) replaces the standard log, yielding a finite running coupling without encountering any Landau pole. In the limit s → 0 the function diverges as S(s) ∼ (1/s) and the coupling goes to zero as a power α(s) ∼ s. A maximum is found at the point where dS(s)/ds = 0, which occurs at s M = 1.044. Of course, this point does not depend on any parameter and provides an independent way to fix the scale m by a comparison with the lattice. From the data of Ref. [24] in Fig.10 the maximum occurs at p ≈ 0.6 − 0.7 GeV yielding a scale m ≈ 0.6 GeV, not too far from the values in Table I. Taking m = 0.6 GeV and the maximum as renormalization point, namely p = 0.67 GeV and α s = 1.21 as in Fig. 13, the plot of Eq.(65) is shown in Fig.10 as a broken line. The running coupling α(s) in Eq.(65) provides a nice qualitative description from first-principles, incorporates the standard leading UV behaviour at one-loop and can be used for extending the standard one-loop running coupling deep in the infrared.

VIII. DISCUSSION
Let us summarize the main findings of the paper. It has been shown that, from first principles, without changing the original Lagrangian, Yang-Mills theory can be studied by a perturbative expansion by just taking a massive propagator as the expansion point. Without the need to include spurious parameters or mass counterterms, the expansion can be renormalized and all the divergences are canceled by the standard wave function renormalization of the fields.
At one-loop, the derivatives of the inverse propagators are determined, up to irrelevant multiplicative factors, by the derivatives of the universal functions F (s), G(s), that do not depend on any parameter. Thus, once a scale is fixed (the theory does not contain a scale that must come from the phenomenology), the inverse dressing functions are determined up to an integration constant. The relevant features of the dressing functions are contained in the universal functions F , G regardless of the specific value of the bare coupling and of N . That scaling property has been shown to be satisfied very well by the lattice data, enforcing the idea that the infrared range of QCD can be studied by perturbation theory.
While the derivatives of the dressing functions are derived exactly, the propagators depend on the integration constants F 0 , G 0 . If the coupling is small and the oneloop approximation is under full control we would expect that a slight change of the additive constants could be compensated by a change of the irrelevant multiplicative factors. Actually, that only occurs in the UV and deep in the infrared where the effective running coupling is small. In the range 0.5 − 1 GeV, where the coupling reaches its maximum, the propagators are sensitive to the choice of the additive constants. That seems to be a sign that higher loops might be relevant when the effective coupling is larger. Thus, we expect that the sensitivity should decrease when higher loops are included in the calculation.
Even where the coupling α s is not very small and twoloop corrections seem to be relevant, the one-loop calculation may acquire a variational meaning. The dependence on the renormalization constants is a consequence of an overall dependence on the ratio of the two energy scales: the mass parameter m and the renormalization point µ. Since the exact result should not depend on that ratio, the dependence is expected to decrease when higher loops are included in the calculation. Thus a best choice for that ratio could be obtained by some station-ary condition on the observables, requiring that the sensitivity should be minimal in the predicted phenomenology. However, there is no proof that a best choice of the renormalization constants does exist, mimimizing twoloop corrections everywhere. Thus, it is encouraging to know that, by tuning the additive constants, the one-loop calculation already provides an excellent description of the lattice data for the propagators and the running coupling. We conclude that, while not anomalously small in general, two-loop corrections can be minimized by a best choice of the constants.
Moreover, the sensitivity to the additive constants F 0 , G 0 seems to be really negligible below 300 MeV and above 1.5 GeV, namely when α s < 0.6. In those ranges the running coupling collapses on the lattice data without the need to tune any constant or parameter, from first principles and by a fully analytical description.
The existence of an energy range, deep in the infrared, where the one-loop approximation seems to be under full control, could open the way for a more general analytical study of QCD below Λ QCD where many interesting phenomena still suffer the lack of a full description from first principles.  [31,32]. The derivatives are straightforward and have been checked by a software package. The result is where the logarithmic functions L x are and the rational parts R x are