The gluon propagator in linear covariant $R_\xi$ gauges

Explicit analytical expressions are derived for the gluon propagator in a generic linear covariant $R_\xi$ gauge, by a screened massive expansion for the exact Faddeev-Popov Lagrangian of pure Yang-Mills theory. At one-loop, if the gauge invariance of the pole structure is enforced, the gluon dressing function is entirely and uniquely determined, without any free parameter or external input. The gluon propagator is found finite in the IR for any $\xi$, with a slight decrease of its limit value when going from the Landau gauge ($\xi=0$) towards the Feynman gauge ($\xi=1$). An excellent agreement is found with the lattice in the range $0<\xi<0.5$ where the data are available.


I. INTRODUCTION
Almost all the visible mass in the universe arises from dynamical mass generation, a mechanism that converts chiral current quarks into constituent quarks, each carrying one third of the proton mass. The mechanism can be understood as the effect of low-energy gluon clouds dressing the current quark, so that the study of the gluon propagator in the IR becomes of paramount importance for a full comprehension of the mass generation [1][2][3][4][5][6][7]. Unfortunately, even in the pure gauge sector, perturbation theory breaks down in the IR and the results of lattice simulations [5,[8][9][10][11][12][13][14][15][16][17][18] are regarded as the only benchmark for the continuum approaches [19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38] that have been developed. Among them, a purely analytical method has been proposed in the last years [39][40][41], which is based on a change of the expansion point of ordinary perturbation theory and provides explicit and very accurate expressions for the gluon propagator in the Landau gauge [42]. The method relies on a screened massive expansion, with massive propagators in the internal gluon lines of Feynman graphs, and is derived from the exact Faddeev-Popov Lagrangian of pure Yang-Mills theory, from first principles, without adding any phenomenological parameter. The expansion can be seen to emerge from the Gaussian effective potential [43,45] which provides a simple argument for the dynamical mass generation of the gluon and has been also studied at finite temperature [44,45].
In this paper, the massive expansion is extended to the more general case of a linear covariant R ξ gauge and explicit analytical expressions are provided for the gluon propagator at any generic value of the gauge-fixing parameter ξ, yielding new insight into the gauge dependence of the propagator, that cannot be extracted by any other method.
Exploring the gauge dependence of the gluon propagator is in itself important in order to individuate the properties that are gauge invariant and might be directly related to physical observables. Despite that, the covariant R ξ gauge, which is under control at the perturbative level, is basically unexplored in the IR because of con-vergence problems in lattice calculations [46][47][48]. Quite recently, a lattice simulation has been extended up to ξ = 0.5 [49], predicting a saturation of the propagator deep in the IR, with very small deviations from the results in the Landau gauge, but in strong disagreement with some recent predictions of a continuum study [25]. On the other hand, the lattice data seem to be in qualitative agreement with the picture emerging by Nielsen identities in Ref. [22]. Out of the Euclidean space, no information has been reported so far about the analytic properties in R ξ gauge.
On general grounds, because of Nielsen identities [50], we know that the poles and the residues of the gluon propagator, i.e. the principal part, must be gauge parameter independent [51][52][53]. While no information on the existence and properties of the poles can be extracted from lattice calculations in the Euclidean space, the massive expansion provides explicit analytical expressions that can be continued to the complex plane [41]. Some attempts at reconstructing the spectral functions from the lattice data have been reported [54,55] and are in qualitative agreement with the predictions of the expansion [41].
At one-loop, by the massive expansion, a pair of complex conjugated poles were found in the Landau gauge [44], as also predicted by different phenomenological models [56][57][58], again in strong disagreement with other continuum studies [59] based on the truncation of an infinite set of Dyson-Schwinger equations. While the genuine nature of the poles was already shown by studying their behavior at finite temperature [44], their explicit gauge invariance would provide further evidence that they are not artifact of the expansion. Strictly speaking, by changing the expansion point, the Becchi-Rouet-Stora-Tyutin (BRST) symmetry of the quadratic part of the Lagrangian is broken in the expansion and we should not expect that the pole structure might be exactly gauge invariant at any finite order. However, since the total Lagrangian is not modified, the gauge parameter independence must be recovered if the expansion provides a very good approximation of the exact propagator. Thus the gauge parameter independence of the pole structure would give a quantitative estimate of the accuracy in the complex plane, where no comparison with the lattice can be made.
By the same argument, the massive expansion can be optimized by enforcing the gauge parameter independence of the whole pole structure, yielding a fully selfcontained calculation from first principles, without any adjustable parameter or external input. Moreover, once optimized in the complex plane, the result is found in excellent agreement with the lattice data in the Euclidean space, not only in the Landau gauge, but for the whole range, up to ξ = 0.5, that has been explored in the lattice so far [49]. No dramatic difference is found for larger values of ξ and even in the Feynman gauge the gluon propagator is finite in the IR, with a slight suppression of its saturation value compared to the Landau gauge. Being gauge parameter independent, the principal part of the propagator might be directly related to physical observables like glueball masses, as recently discussed by a quite general method [58].
The paper is organized as follows: in Section II the massive expansion of Refs. [39,40] is extended to a generic R ξ gauge; in Section III the expansion is optimized by requiring that the pole structure is gauge parameter independent as demanded by Nielsen identities; in Section IV the optimized gluon propagator is shown for a wide range of the gauge parameter ξ, including the Feynman gauge (ξ = 1), and is compared with the available lattice data; Section V contains a brief discussion of the main results. Explicit analytical expressions for the propagator in R ξ gauge are derived in Appendix A with many details on the calculation of the graphs.

II. THE MASSIVE EXPANSION IN R ξ GAUGE
The massive expansion has been first developed in Refs. [39,40] and related to the Gaussian effective potential in Refs. [43,45]. It is based on a change of the expansion point of ordinary perturbation theory for the exact gauge-fixed Faddeev-Popov Lagrangian of pure Yang-Mills SU (N ) theory. The Lagrangian can be written as where L Y M is the Yang-Mills term the tensor operatorF µν iŝ L F P is the ghost term arising from the Faddeev-Popov determinant and L f ix is the covariant gauge-fixing term The gauge field operators arê where the generators of SU (N ) satisfy the algebra In the standard perturbation theory, the total action is splitted as S tot = S 0 + S I where the quadratic part can be written as and the interaction is with the three local interaction terms that read In Eq. (7), ∆ 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 where η µν is the metric tensor. The massive expansion is obtained by adding a transverse mass term to the quadratic part of the action and subtracting it again from the interaction, leaving the total action unchanged.
In some detail, 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. In order to leave the total action unaffected by the change, we must add the same term in the interaction, providing a new interaction vertex δΓ. Dropping all color indices in the diagonal matrices and inserting Eq.(10) and (14) in Eq.(13) the vertex is just the transverse mass shift of the quadratic part δΓ µν (p) = m 2 t µν (p) (15) and must be added to the standard set of vertices in Eq. (9). The proper gluon polarization Π and ghost self energy Σ can be evaluated, order by order, by perturbation theory. In all Feynman graphs the internal gluon lines are replaced by the massive free-particle propagator ∆ m µν and all insertions are considered of the (transverse) mass counterterm δΓ µν which plays the role of a new two-point vertex. It is shown as a cross in Fig.1 where some twopoint self-energy graphs are displayed. We will refer to the graphs with a cross as crossed graphs.
Since the total gauge-fixed FP Lagrangian is not modified and because of gauge invariance, the longitudinal polarization is known exactly and is zero, so that the total polarization is transverse and the (exact) dressed propagators read where the transverse and longitudinal parts are At tree level, the polarization is just given by the counterterm δΓ of Eq. (15), so that the tree-term Π tree = m 2 just cancels the mass in the dressed propagator ∆ of Eq.(18), giving back the standard free-particle propagator of Eq.(10).
Finally, summing up the loops and switching to Euclidean space, the transverse dressed propagator can be written as where Π loops (p) is given by the transverse part of all the loop graphs for the (proper) polarization.
At one-loop, as discussed in Refs. [39,40], we sum all the graphs with no more than three vertices and no more than one loop, which are displayed in Fig. 1. In Appendix A, explicit analytical expressions are given for all the polarization graphs of the figure.
The diverging integrals are made finite by dimensional regularization and can be evaluated in the Euclidean space, by setting d = 4 − ǫ. An important feature of the massive expansion is that the crossed graphs cancel all the spurious diverging mass terms exactly, so that no mass renormalization is required. That is a very welcome feature since there is no bare mass in the original Lagrangian. At one-loop, as shown in Appendix A, in the M S scheme, the diverging part of the proper transverse polarization can be written as which is the same identical result of standard perturbation theory [60] and ensures that we obtain the correct leading behavior in the UV where the mass insertions are negligible, as shown in Eq.(A49). As usual the diverging part can be canceled by wave function renormalization, by subtraction at an arbitrary point. Of course, a finite term ∼ const. × p 2 arises from the subtraction and cannot be determined in any way. It also depends on the regularization scheme and on the arbitrary scale µ, so that its actual value remains somehow arbitrary. It basically is the only free parameter of the approximation, as discussed later. For an observable particle, the constant would be fixed on mass shell, by requiring that the pole of the propagator is at the physical mass with a residue equal to 1. The confinement of the gluon has been related to the existence of complex conjugated poles [44], so that if, on the one hand, there is nothing like an observable gluon mass, on the other hand, the analytic properties at the poles and their gauge parameter independence will be shown to be enough for determining the propagator entirely and uniquely.
The finite part of the one-loop proper polarization, as resulting from the sum of all the graphs in Fig. 1, reads Figure 1: Two-point self-energy graphs with no more than three vertices and no more than one loop.
where s = p 2 /m 2 is the Euclidean momentum. The functions F (s) and F ξ (s) are adimensional and do not depend on any parameter. Their explicit expressions are derived in Appendix A by a detailed calculation of the integrals and the final result is reported in Eqs.(A41),(A44). The constant C arises from the subtraction of the diverging part by wave function renormalization. For a generic subtraction point p = µ, the one-loop transverse propagator follows from Eq.(19) where Z µ is the arbitrary finite renormalization constant Z µ = µ 2 ∆(µ). Finally, the propagator can be written as where the coupling and all other constants are absorbed by a finite renormalization factor Z and the new constant F 0 which depend on the subtraction point µ according to Eq.(23) provides an explicit analytical expression for the one-loop gluon propagator. It contains three parameters: m, Z and F 0 . However, the finite renormalization factor Z is irrelevant, while m is the unique energy scale. Since the exact Lagrangian does not contain any energy scale, m cannot be determined by the theory: the mass parameter m determines the overall energy scale and can only be fixed by comparison with some physical observable. That is not a limitation of the approximation but is a standard feature of Yang-Mills theory. Moreover, being just a scale parameter, the mass m is not a physical or dynamical mass and is not even required to be gauge invariant. We will use the energy scale of the lattice and fix m by comparison with the data of simulations in the Landau gauge. Thus, the only free parameter in Eq. (23) is the constant F 0 which is related to the arbitrary ratio µ/m. Since the result does depend on F 0 , the expansion must be optimized by a criterion for determining the best F 0 , yielding a special case of optimized perturbation theory by variation of the renormalization scheme, a method that has been proven to be very effective for the convergence of the expansion [61].
Assuming that the expansion converges more quickly for an optimal value of F 0 , the one-loop result might be very close to the exact result for a special choice of the constant. That is shown to be the case in Refs. [39][40][41][42] where an excellent agreement with the lattice is found in the Landau gauge. Unfortunately, the available data are not fully consistent and a best fit yields slightly different values of F 0 and m for different data sets, as shown in Table I. The deviations might be related to a slightly different choice of units as recently discussed in Ref. [62]. We can extract a global average F 0 ≈ −0.9 ± 0.1. Of course, the actual value of the constant F 0 depends on the details of the definition of the functions F (s), F ξ (s) which are evaluated up to an (omitted) arbitrary additive constant in Appendix A. In this paper, all the values of F 0 refer to the definition given by Eqs.(A41),(A44) for those functions.   [18] (in the range 0 -4 GeV) and Ref. [12] (0 -2 GeV), and by the SU (2) data of Refs. [9,10] (0 -2 GeV).
In the Landau gauge, the best agreement is found for the data set of Ref. [18], with a best fit parameter F 0 = −0.887 and a mass scale m = 0.654 GeV, evaluated in the range 0 < p < 4 GeV. The resulting gluon propagator is shown in Fig. 2 together with the lattice data.

III. OPTIMIZATION FROM FIRST PRINCIPLES
If the expansion is optimized in the Euclidean space, by a direct comparison with the lattice data, any control of the approximation is lost in Minkowski space and one might wonder how robust the optimal choice would be when continued to the complex plane. Moreover, a selfcontained optimization strategy, which does not require any external input, would be essential for exploring new aspects that are out of the reach of lattice calculations. In this section, we show that the expansion can be optimized from first principles in the complex plane by enforcing some general exact analytic properties that arise from the BRST invariance of the gauge-fixed Lagrangian.
The Nielsen identities [50] are exact equations connecting the gauge parameter dependence of some correlation functions with other Green functions. Their proof follows from the BRST invariance of the Faddeev-Popov Lagrangian, Eq.(1), which has not been modified by our change of the expansion point. They have been used as a tool for establishing general invariance properties of the pole structure in QCD [52] and in other Yang-Mills theories [57].
Following the detailed derivation of Ref. [52], the exact transverse projection of the gluon propagator ∆(p) must satisfy the Nielsen identity where, omitting the diagonal color indices, G T (p) is the transverse component of the Green function G µν ab (−p, p, 0) which is defined as in terms of the Nakanishi-Lautrup auxiliary field B a and of the covariant derivative of the ghost field D µ ω a . If the gluon propagator has a pole in the complex plane at p 2 = p 2 0 (ξ), then the inverse propagator has a zero and we can write the identities Then, the vanishing of the right hand side of Eq.(25) at p = p 0 (ξ) says that the partial derivative is also zero and the pole p 0 must be gauge parameter independent By the same argument, the residues at the poles are also gauge parameter independent [53]. In fact, if we differentiate Eq.(25) with respect to p 2 ∂ ∂ξ d dp 2 the right hand side vanishes at p = p 0 because of Eq. (28), so that the residue R, defined as satisfies the exact equation We conclude that, for the gauge-fixed Yang-Mills Lagrangian, the principal part ∆ P of the exact gluon propagator must be gauge parameter independent. The argument fails if G T (p) has a pole in p = p 0 , which is usually not the case.
In the quadratic part of the Lagrangian, the BRST symmetry is broken by the mass term that has been added and has been subtracted again from the interaction. Thus, while the total Lagrangian is BRST invariant, the symmetry is broken at any finite order of the massive expansion. For that reason, we do not expect that the one-loop propagator might satisfy the Nielsen identity exactly. However, the closer we reach to the exact result, the better is expected to be the agreement with the exact identities. Thus, we can exploit the dependence on the parameters F 0 , m in Eq.(23) and optimize the expansion by requiring that the pole structure of the propagator is gauge parameter independent. That is equivalent to an optimal choice of the subtraction point µ/m, which is usually fixed on mass shell for an observable particle. Without any observable gluon mass at hand, the invariance of the poles and residues turns out to be enough for determining the one-loop gluon propagator entirely and for any choice of the gauge parameter.
For a generic choice of the gauge parameter ξ, the optimal parameters can be regarded as functions F 0 (ξ), m(ξ), to be determined by the requirement that the pole and the residue do not depend on ξ. Of course, the finite renormalization factor Z remains arbitrary and has no physical relevance. Let us denote by Ψ(z, ξ, F 0 , m) the inverse dressing function in Eq.(23) which is an analytic function of the complex variable z = x + iy. is the Euclidean momentum. On the real axis, for y = 0, we recover the Minkowskian momentum p 2 M = z 2 = x 2 . Thus, the variable z is the analytic continuation of the physical momentum p M . The pole z 2 0 = −p 2 0 is a zero of the inverse dressing function Ψ and must satisfy the equation Ψ(z 0 , ξ, F 0 , m) = 0. The gauge parameter independence of the pole requires that yielding a set of two coupled real equations for the real and imaginary parts. The equations can be solved for F 0 (ξ 2 ) and m(ξ 2 ) from a given initial value F 0 (ξ 1 ), m(ξ 1 ). Taking the Landau gauge as the initial point ξ 1 = 0 and fixing a scale m 0 = m(0) as energy units, the functions F 0 (ξ) and m(ξ) are determined for any value of the gauge parameter ξ from the initial value F 0 (0) which remains the only free parameter. Thus, we can encode the gauge parameter independence of the pole in the optimized propagator and evaluate it for any value of the parameter ξ. The functions F 0 (ξ), m 2 (ξ) are shown in Fig. 3 and Fig. 4 for different choices of the initial value F 0 (0) in the Landau gauge.
In the range −2 < F 0 (0) < 0, the gluon propagator of Eq.(23) has a single pair of complex conjugated poles, while other values of F 0 (0), out of that range, seem to be unphysical. For F 0 (0) < −2 the expression in Eq. (23) has poles in the Euclidean space and changes sign at the poles, on the positive s axis. Moreover, according to Eq.(24), the coupling g 2 would become negative in that range because the minimal value of F (s) is ≈ 2. For F 0 (0) > 0 the coupling g 2 becomes very small in Eq. (24) for any µ and the pole topology becomes very different. As discussed in the previous section, in the Landau gauge, the best agreement with the lattice is found for F 0 ≈ −0.9 which is at the center of the physical allowed range.
It is remarkable that, close to the best fit value F 0 (0) ≈  Table  II). The curves are approximately tangent (i.e. θ ≈ 0) at the intersection point z0 (the pole) whenever F0 ≈ −0.9.
lines in the map and setting ξ 1 = 0, ξ 2 = ξ, we can write and because of Eq.(31), the angle θ gives the phase change of the residue R which can be written, as a function of ξ, since the modulus |R| can always be made invariant by an appropriate choice of the real renormalization constant Z(ξ). Explicit analytical expressions for the derivative of Ψ are reported in Eqs.(A50),(A53) of Appendix A. We observe that the angle θ is not exactly zero, so that in general, the Nielsen identity Eq.(25) and its consequences Eqs. (29),(32) cannot be all satisfied. However, as shown in Fig. 5 and Fig. 6, the angle θ becomes very small, for a wide range of ξ, if the initial constant F 0 (0) is close to the value F 0 ≈ −0.9 which already described the lattice data very well in the Euclidean space. In other words, the optimal propagator in the Euclidean space is also the one that best satisfies the Nielsen identity in the complex plane, giving us confidence in the general accuracy of the approximation. We must mention that averaging over Gribov copies might break BRST invariance in the lattice. However, we are assuming that the Nielsen identities are not seriously affected in lattice calculations.
Reversing the argument, the expansion can be optimized in a self-contained way, by first principles and without any external input, by assuming that the best choice for the initial constant F 0 (0) is the one that makes the angle θ smaller in a wider range of ξ. Even if there are no technical reasons for limiting the value of the gauge parameter, we expect that perturbation theory would be more effective when ξ is small and the expansion might be out of control for very large ξ ≫ 1. Prudentially, the present study is limited to the range ξ < 1.2, including the Feynman gauge.
The minimal phase deviation is observed for the initial value F 0 (0) = −0.876. As shown in Fig. 6, for that choice, the phase θ fluctuates around zero in the whole range 0 < ξ < 1.2, with very small deviations which are less than 0.003. Nevertheless, no fine tuning is required since θ is very small around F 0 (0) ≈ −0.9 and any slight change of F 0 (0) can be compensated by an appropriate choice of Z and m. In fact, as shown in Fig. 2, when the present new set of first-principle optimal parameters are inserted in Eq.(23), the propagator is indistinguishable from the previous one that was obtained by a best fit of the lattice data. Actually, the optimal initial value F 0 (0) not only minimizes the phase deviation θ(ξ), but also makes m 2 (ξ) stationary and maximal for any fixed ξ, as shown in Fig. 3 where the optimal curve is plotted as a red line. That is a geometric consequence of the pole being the tangency point in Fig. 5. Moreover, the optimal function F 0 (ξ) is the most gauge parameter invariant curve in Fig. 4 (shown as a red line).
The optimal parameters are summarized in Table II together with very accurate polynomial interpolation formula for the optimal functions F 0 (ξ), m 2 (ξ). Extracting the energy scale m 0 = m(0) = 0.656 GeV from the lattice data of Ref. [18] in the Landau gauge, the invariant pole is found at x 0 = M = 0.581 GeV and y 0 = γ = 0.375 GeV, which might be regarded as the physical mass and the damping rate of the quasigluon, respectively, as discussed in Ref. [44].

IV. THE PROPAGATOR AT ξ = 0
In the Euclidean space, the gluon propagator can be evaluated analytically by Eq.(23), for any value of the gauge parameter ξ, inserting the optimal parameters of Table II which Table II: Set of optimal parameters, obtained by enforcing the gauge parameter independence of the pole structure in the range 0 < ξ < 1.2. The energy scale m0 and the finite renormalization constant Z(0) are determined by the data of Ref. [18] which are shown in Fig. 2.  Table II, for ξ = 0, 0.5, 1 and renormalized at µ = 4.317 GeV. The points are the lattice data of Ref. [49].
compare with the available lattice data of Ref. [49], the finite renormalization constant Z is fixed by the same momentum subtraction scheme of that work, i.e. requiring that µ 2 ∆(µ) = 1 for any ξ and taking the same renormalization point µ = 4.317 GeV. That is equivalent to taking the constant Z in Eq.(23) to be Z = Ψ(iµ, ξ, F 0 , m).
The gluon propagator is shown in Fig. 7, for several values of the gauge parameter ξ, together with some data points extracted from Ref. [49].The agreement with the data is very good in the limited range ξ < 0.5 where they are available. For ξ = 0, the propagator is slightly suppressed in the IR compared with the Landau gauge. We must mention that previous continuum studies, based on the truncation of an infinite set of exact Dyson-Schwinger equations, reached contrasting and ambiguous results. While a strong dependence on the gauge parameter was predicted in Ref. [25], with large deviations from the Landau gauge, a qualitative agreement with the lattice was  Table II, for ξ = 0.1, 0.5 and 1. The bars are the lattice data of Ref. [49].
reported in Ref. [22] by the aid of exact Nielsen identities which seem to play a key role. The gauge dependence was found small but no quantitative prediction could be made and even the sign of the change was not defined by that method. As shown in Fig. 7, up to and beyond the Feynman gauge (ξ = 1), no dramatic change occurs and the suppression of the propagator increases very smoothly with the increasing of ξ. The change can best be seen by evaluating the ratio between ∆(p) at ξ = 0 and at ξ = 0, as shown in Fig. 8 together with the lattice data of Ref. [49]. Even if the lattice calculation is plagued by large statistical errors, with scattered data and large error bars, the optimized propagator seems to be in quantitative agreement with the data and reproduces the correct trend predicted by the lattice. We stress that the curves are not a fit of the data and the agreement is reached from first principles without any adjustable parameter. The dressing function is shown in Fig. 9. As predicted by the lattice [49], the maximum is basically fixed at the same energy for any ξ. We argue that the Nielsen identity gives the correct scale factor m(ξ)/m(0) that keeps the maximum fixed, at variance and in strong contrast with the continuum calculation of Ref. [25] which might miss that important constraint.
In the studied range of ξ, the whole principal part of the propagator in Eq.(33) is basically invariant up to a finite renormalization factor. The pole p 0 is fixed at the value of Table II, while the phase of the residue is Arg R(ξ) = 1.262 + θ(ξ) where |θ(ξ)| < 2.75 · 10 −3 , yielding the ratio t R = Im R(ξ)/ Re R(ξ) = 3.132 ± 0.03. This ratio is important for determining the explicit parameters of the rational part Eq.(33) which has been derived at tree level by other phenomenological models like the refined Gribov-Zwanziger model [56][57][58]. Being gauge parameter independent, the parameters of the rational part might be directly related to physical observables or condensates [63,64] and a recent general method has been proposed for extracting information on the glueball masses [58]. Using the notation of Ref. [64], the principal part of the propagator, Eq.(33), can be written as where having made use of the optimized parameters of Table  II. Below 1 GeV, the masses M i seem to be compatible with the statistical analysis of Ref. [64], even if the simple rational part ∆ P was used in that work for a fit of the lattice data, ignoring the corrections which are included in the present optimized one-loop propagator. In fact, the corrections are gauge dependent and very small below 1 GeV, as already shown in the Landau gauge by a direct evaluation of the spectral function [41,65]. The Schwinger function ∆(t) can be evaluated by a numerical integration, as a function of the Euclidean time t, according to its definition and is shown in Fig. 10 for different values of the gauge parameter. In the Landau gauge, the Schwinger function is found in qualitative agreement with the result of Ref. [66], with a positivity violation that occurs above the point t = t 0 ≈ 5.8 GeV −1 where the function crosses the zero and becomes negative. The scale t 0 is roughly the size of a hadron and in Ref. [66] it was conjuctered to be a physical gauge-invariant scale at which gluon screening occurs. Actually, as shown in Fig. 10, the crossing point t 0 is found to be almost gauge parameter independent. Moreover, the large t behavior seems to be dominated by the singularities and the whole Schwinger function is very well approximated by inserting in Eq.(41) the simple principal part ∆ P (p) of Eq. (33), which is gauge parameter independent, yielding the analytical result ∆ P (t) = |R| M 2 + γ 2 e −Mt cos γt − θ + arctan γ M (42) which is shown in Fig. 10 as a broken line.
We cannot end this section without a brief discussion of the spectral function, which has attracted great interest [54,55] even if its physical content is quite unclear in presence of complex poles and confinement. In fact, the usual Källen-Lehmann representation must be replaced by the more general integral representation [65] Re ∆(p 2 ) = ∆ P (p 2 ) + P.V.
where the spectral function ρ(p 2 ) is gauge dependent and does not contain any information on the gauge parameter independent principal part ∆ P which must be added to the integral for reproducing the whole propagator. Moreover, ρ(p 2 ) is even not positive defined for a confined particle. In the Landau gauge, the spectral function was evaluated by the massive expansion in Ref. [41] and the dispersion relation of Eq. (43) was checked in Ref. [65] by a numerical integration. The integral provides the difference between the principal part and the whole propagator, so that the difference can be large only if the total weight which comes from the integration of ρ(p 2 ) is large. Moreover, ρ(p 2 ) changes sign and the contributions arising from different signs can partially cancel. The one-loop spectral density can be easily evaluated by the explicit expression of the propagator, Eq.(23), using the optimal parameters of Table II, and is shown in Fig. 11 for different values of the gauge parameter ξ. It has some gauge dependent features, like a cusp at the two-particle threshold p = 2m(ξ) and a finite spike at p ≈ m(ξ). In the Landau gauge, the spike is just a smooth maximum but is enhanced for ξ > 0.08 by the appearance of a gauge dependent pole near the real axis, at x ≈ m(ξ). Some details of the finite peak on the real axis are shown in Fig. 12. Apart from the peak, the spectral density is very small and even the peak area gives a small contribution to the integral in Eq.(43) because of the change of sign that occurs just at the peak, in agreement with a confinement scenario. While the peak resembles the spike which was observed in Ref. [59], its physical nature is unclear and is certainly related to the nature of the new gauge dependent pole which might be an artifact of the one-loop approximation.
From a technical point of view, the new pole arises because of the logarithmic divergence of the real part of F ξ (−z 2 /m 2 ) at the branch point x = m on the real axis. The divergence occurs because of the bad IR behavior of the crossed gluon loop, Π 2c in Fig. 1, in the limit p → im, since the denominator in Eq.(A19) becomes point, for x ≈ m and any finite ξ = 0, the real part of the inverse dressing function Ψ(z), Eq. (34), can be written as where Ψ reg is the regular part and the prefactor A(z) of the log is a rational function which is real on the real axis, with A(m) = −2/3. Then, taking z = m + re iφ , the contour line Re Ψ = 0 is given by which is a very small circle centered at x = m on the real axis, with an exponentially small radius in the limit ξ → 0 if C > 0. In the Feynman gauge, ξ = 1, the contour line is just visible in Fig. 5 as a small black semicircle centered at x = m(1) = 0.53 GeV on the real axis. It gets hardly visible for ξ < 0.5. At the same branch point, the imaginary part of F ξ (s) has a large discontinuous step yielding a change of the whole imaginary part which is quite larger than Im Ψ reg (m) ≈ 0.17 and gives rise to a sharp change of sign at x = m(ξ), even when ξ is small, provided that ξ > 0.08. On the complex plane, the discontinuous step is smeared out and the imaginary part Im Ψ changes sign on a contour line Im Ψ = 0 which originates from the branch point x = m(ξ), just at the center of the circle Re Ψ = 0. The resulting contour line Im Ψ = 0 is visible in Fig. 5 as a solid red line ending at the center of the black semi-circle. The crossing point of the two contour lines is the new pole that appears for ξ > 0.08. On the other hand, if ξ < 0.08, the imaginary part Im Ψ changes sign below x = m(ξ), out of the circle, the contour lines do not cross and the extra pole disappears when approaching the Landau gauge. By the previous analysis we conclude that the narrow peak of the spectral function must have a very small width, roughly given by the distance of the pole from the real axis r ≈ exp(−C/ξ), getting smaller and smaller when ξ ≪ 1, as shown in Fig. 12. Moreover, Im Ψ and ρ(p 2 ) change sign across the peak and the overall effect of the peak on the integral, in Eq. (43), is expected to be negligible.
It is likely that the sharp peak of the spectral function and the gauge dependent pole get smoothed in the exact propagator since the Nielsen identity, Eq.(25) would forbid the existence of a pole which depends on the gauge parameter ξ, unless the Green function G T (p) in Eq. (26) has a pole at the same point. Having traced the source of the pole and found it related to the logarithmic divergence of the crossed graphs at p 2 = −m 2 , we cannot exclude that the same divergence might occur in the ghost sector and in other Green functions. Thus, in principle, we cannot rule out that the pole might be genuine, even if probably related to unphysical degrees of freedom of the ghost sector.

V. DISCUSSION
There is a growing consensus that QCD and Yang-Mills theory are self-contained theories that dynamically generate their own infrared cutoff. The numerical simulations on the lattice have shown that the exact theory generates a dynamical mass which screens the gluon interaction in the IR. Therefore, any continuum first-principle study should reproduce the same results without the aid of any adjustable parameter, except for the overall energy scale that must come from the phenomenology. It could be argued that, because of Gribov ambiguity, in R ξ gauge the Faddeev-Popov Lagrangian is just an approximation of the full theory. The approximation works very well in the usual perturbative approach but could be out of control in the IR because of non-perturbative effects. A phenomenological parameter has been introduced by several authors for locating the Gribov horizon, yielding an interaction-induced mass scale which screens the theory in the IR [56][57][58][67][68][69][70][71][72]. However, even averaging over Gribov copies, a dynamical mass is generated in the theory, as shown by the gauge-fixed lattice calculations in the Landau gauge. A recent analysis [73] has made clear that the dynamical mass would be as effective as the Gribov parameter for screening the theory and that its dynamical appearance alone would eliminate the problem of Gribov copies and complete the definition of the theory.
The same argument holds for the massive expansion which is a screened expansion from the beginning and can be safely used in the IR. Having changed the expansion point, the gauge-fixed theory can be studied by plain perturbation theory and the agreement with the lattice data shows that, when the expansion is optimized, higher order graphs are very small and negligible. Thus, ignoring the Gribov ambiguity does not seem to be a problem as far as perturbation theory works well. Again, it is a consequence of the dynamical mass that screens the theory, yielding a self-contained perturbative description from first principles.
It is not surprising that, without using any adjustable parameter and without modifying the original gaugefixed Lagrangian, the massive expansion predicts the same pole structure which was found by the refined Gribov-Zwanziger model [56][57][58]. The two approaches are very different but they study the same identical physical system, so that if both are valid approximations they must reach the same conclusions. Moreover, our analysis supports the physical relevance of the principal part: having established its gauge parameter independence [53], we argue that the simple rational part ∆ P (p) might play an important role in the phenomenology, more than the (small) gauge dependent spectral density. The conclusions of the present work would be enforced by a comparison with position-space lattice data, because of their sensitivity to the analytical structure of the propagator. Unfortunately, at the moment, for a generic covariant gauge, no such data are available.
An apparent drawback of the massive expansion is that the BRST invariant action is arbitrarily splitted in two parts that are not BRST invariant. The Nielsen identities cannot be satisfied exactly at any finite order of the expansion. However, because of the spurious dependence of the approximation on the subtraction point µ/m, the expansion can be optimized by enforcing the gauge parameter invariance of the pole structure. Thus, the extension to R ξ gauge, not only gives new information on the gluon propagator in a generic gauge, but also provides a unique way to fix the optimal expansion even in the Landau gauge. The good agreement with the available lattice data, which is reached without any fit of adjustable parameters, increases our confidence in the general validity of the method as a first-principle benchmark for more phenomenological models. vertex which is shown as a cross in the figure. We refer to the graphs with one insertion of the counterterm as crossed graphs.
1. Graphs Π 1b , Π1c and Π 1d (tadpoles) In the Euclidean space, the constant tadpole Π 1b can be written as having dropped the longitudinal loop which is scaleless and vanishes in dimensional regularization. Setting d = 4 − ǫ, in the M S scheme, where C is a constant which depends on the regularization scheme. The crossed graphs do not contain any longitudinal gluon line since the counterterm δΓ is transverse in Eq. (15). The graph Π 1c can be written as a derivative (A3) As expected, the diverging terms cancel in the sum Π 1b + Π 1c . The double-crossed tadpole Π 1d is finite and including its symmetry factor it reads so that the sum of the constant graphs is The ghost loop Π 2a is a standard graph and does not depend on ξ. In the Euclidean space it is given by the integral [32] The integral is straightforward and setting d = 4 − ǫ the diverging part is while the finite part reads where s = p 2 /m 2 and the constant C 0 depends on the regularization scheme.

Gluon loop Π 2b
The gluon loop Π 2b can be written as where Π 0 2b (p) is the graph in the Landau gauge, ξ = 0. In the Euclidean space, setting d = 4, it reads [32] Π 0 2b (p) = N g 2 6 and the kernel F 0 can be derived by the explicit expressions of Ref. [32] F 0 (k, p) = 10(k 2 + p 2 ) + (k + p) 2 (A11) It is useful to decompose it as and using the identity the graph can be split as where The integrals can be evaluated analytically [40,67,68] by dimensional regularization for d = 4 − ǫ, yielding a diverging part and a finite part where C 1 , C 2 are constants which depend on the regularization scheme, s = p 2 /m 2 and L A , L B are the logarithmic functions L A (s) = (s 2 − 20s + 12) 4 + s s The other terms, Π ξ 2b and Π ξξ 2b , arise by substituting one and two transverse lines, respectively, with the longitudinal ones. By the general scheme of Ref. [32], for d = 4, they follow as The quadratic term is trivial since the integral Π ξξ 2b is scaleless and by a dimensional argument Π ξξ 2b (p) = const × p 2 . The constant can be absorbed by a finite wave function renormalization and the term can be ignored.
The two integrals in Eq.(A19) must be the same, as can be easily seen by substituting k → (−k − p) in Eq.(A21). Taking twice the explicit expression of F ξ0 , the integral can be written as By dimensional regularization, taking d = 4 − ǫ, the integrals can be evaluated analytically in the M S scheme. The first integral is the same occurring in Eq.(A2) The other integrals are where all constants C X , C ′ X depend on the regularization scheme. In Eq.(A27), the first two lines arise from the p 4 term of I ξ D (p) and the diverging term cancels the corresponding divergence of I ξ C (p) in Eq.(A26). Adding up the different integrals we obtain a diverging part and a finite part where we have omitted the irrelevant constants.
Finally, the gluon loop has the following structure

Standard one-loop graphs
The standard one-loop result of perturbation theory does not contain any contribution from the crossed graphs. In a generic linear covariant gauge, the standard one-loop polarization Π 1 (p) is obtained as the sum and summing up the explicit expressions reported above, we find a diverging part In the limit m → 0 the diverging part in Eq.(A32) agrees with the well known result of perturbation theory [60].
In the limit ξ → 0 the finite part in Eq.(A33) gives the known result in the Landau gauge [67,68]. The constants C m and C p are arbitrary since they depend on the regularization scheme and on the arbitrary energy scale µ in Eq.(A32). In the standard perturbation theory, they are the finite parts resulting from the cancellation of the divergences by mass and wave function renormalization, respectively. In pure Yang-Mills theory, there is no mass term in the original Lagrangian and no mass renormalization for the cancellation. However, all constant mass terms cancel exactly by inclusion of the crossed graphs.

Total polarization (including the crossed graphs)
All crossed graphs, containing one insertion of the transverse mass counterterm, can be added to the total one-loop polarization by a simple derivative, as discussed above, for the tadpole. The sum of all graphs in Fig. 1 follows as where f ′ (s) and f ′ ξ (s) are the derivatives of the functions f (s) and f ξ (s), respectively.
Finally, inserting the polarization in Eq. (19) and canceling the divergence by the usual wave function renormalization, the renormalized dressed propagator reads where Z is an arbitrary finite renormalization factor, F 0 is a finite additive constant and the adimensional functions F , F ξ do not depend on any parameter and are defined as Their explicit expressions follow by the simple derivative of Eq.(A34). The function F (x) was first derived in Refs. [39,40] and it reads where the logarithmic functions L x are and the rational parts R x are R a (x) = − 4 + x x (x 2 − 20x + 12) The explicit expression of F ξ (x) is and has the leading behavior in the limit x → 0 In the same IR limit, the transverse propagator is finite and the mass parameter M 2 ξ is defined as In the limit x → ∞, the asymptotic UV behavior is The discussion on gauge invariance requires the derivatives of the functions F (x) and F ξ (x). The derivative of F (x) reads where the logarithmic functions L ′ x , for x = a, b, c, are c (x) = −6x log x (A51) and R(x) is the sum of all the rational terms coming out from the derivatives The derivative of F ξ (x) reads