Thermal extension of the screened massive expansion in the Landau gauge

The massive screened expansion for pure SU(3) Yang-Mills theory is extended to finite temperature in the Landau gauge. All thermal integrals are evaluated analytically up to an external one-dimensional integration, yielding explicit integral representations of analytic functions which can be continued to the whole complex plane. The propagators are explored in the Euclidean and Minkowski space from first principles and are compared with the available lattice data. The dispersion relations of the quasi-gluon are extracted from the pole trajectory in the complex plane. A crossover is found for the mass, suppressed by temperature like an order parameter in the confined phase, while increasing like an ordinary thermal mass in the deconfined phase. Despite the high accuracy which is reached by the method at zero temperature, deep in the IR the agreement with the lattice data is only qualitative at finite temperature. The description improves by a temperature-dependent tuning of the optimal renormalization scale and of the mass parameter.


I. INTRODUCTION
In the last decades, considerable efforts have been devoted to the study of the complex behavior of quarks and gluons under the extreme conditions which are reached in heavy-ion collisions. In principle, the dynamical and thermal properties of a quark-gluon plasma should descend from the relatively simple Lagrangian of the SU(3) gauge theory which describes QCD. However, things are not so easy because the standard perturbative approach breaks down in the strong-coupling IR limit and is also plagued by further resummation problems at any finite temperature. As a matter of fact, we still miss a full theoretical treatment of the problem.
Even the pure gauge theory, without quarks, is not fully understood, despite its relevance for describing the quark-gluon plasma. Many important advances have been made by the numerical simulation of the pure Yang-Mills (YM) Lagrangian on a lattice, providing insights into the gluon dynamics and the phase diagram. Among them, the confirmation of a dynamically generated gluon mass [1][2][3][4][5][6][7][8], as predicted by Cornwall in 1982 [9], and the occurrence of a phase transition, with the gluons that become confining below a critical temperature [10][11][12].
It would be a desirable progress if the dynamical and transport parameters, like masses, widths, dispersion relations, transport coefficients, etc., which are currently regarded as phenomenological parameters [13][14][15][16], could be directly evaluated from first principles. That program might be accomplished in part if the elementary correlators and their analytic properties were known in the Minkowski space. Unfortunately, all lattice calculations and most numerical works provide information in the Euclidean space and the analytic continuation is a difficult ill-defined problem for the numerical data [17].
In the last years, a very predictive analytical method has been developed [18][19][20][21] by a mere change of the expansion point of ordinary perturbation theory (PT) for the exact gauge-fixed Becchi-Rouet-Stora-Tyutin (BRST) invariant YM Lagrangian, yielding a screened massive expansion which is safe in the IR while recovering the correct results of ordinary PT in the UV. At oneloop and zero temperature, the screened expansion provides analytical results which are in excellent agreement with the lattice and can be easily continued to Minkowski space [21][22][23][24][25]. Thus the method provides a way to extract dynamical details like masses and damping rates from first principles.
In this paper, the formalism is extended to a finite temperature T = 0, with the aim to provide a complementary tool for the study of the gluon plasma from first principles. As briefly discussed by one of us [26], the screened expansion can be extended to finite temperature, providing a quasi-particle picture for the gluon which is damped, with a very short finite lifetime, and canceled from the asymptotic states. Here, we give a full account of the details of the calculation and report a comprehensive set of results for the gluon sector, including propagators, analytic properties, poles, masses, widths and dispersion relations. We discuss different optimization strategies and, by a comparison with the available lattice data, we explore how robust the screened expansion is when it is extended to finite temperature.
While the existence of a screening mass mitigates the effects of the hard thermal loops, several problems arise at a finite temperature, ranging from the temperature dependence of the optimal mass scale, to the analytic continuation of the numerical integrals. Actually, even if a formal extension to finite temperature is straightforward and based on standard thermal Feynman graphs, the ambition to extract analytical results requires a quite tedious and lengthy analytical calculation of the integrals and, even so, a final one-dimensional numerical integration cannot be avoided. Nonetheless, the resulting numerical integrals are shown to define analytic functions which can be evaluated in the complex plane. Then, the poles of the gluon propagator and the resulting dispersion relations can be easily extracted numerically.
Overall, despite the expected difficulties, the one-loop screened expansion seems to be reliable at low temperature, with correct predictions which become less quantitative at high temperature, especially for the longitudinal sector, when compared with the lattice data.
At T = 0, the one-loop approximation is quite sensitive to the renormalization scheme and to the subtraction point, but it can be shown to be basically tangent to the exact result, which is approached for a special choice of the ratio between the gluon mass parameter m and the renormalization scale µ. Here, m is just a mass parameter which defines the shift of the expansion point [18,19,24,25], not to be confused with the physical mass of the gluon. It seems that, for that special ratio µ/m, the higher order terms become negligible, yielding very accurate analytical expressions for the propagators. While that special ratio is scheme-dependent, it can be determined from first principles by monitoring some identities which must be fulfilled by the exact propagators, like the Nielsen identities, which express the gaugeinvariance of the poles [21]. We must mention that, once the ratio is optimized in the complex Minkowski space, where the poles are defined, the propagators are found in excellent agreement with the lattice data in the Euclidean space. Thus, the optimized analytical expression is not just a good interpolation formula, but a very good approximation for the whole analytic function which is defined in the complex plane. Moreover, at the optimal ratio µ/m there is only one energy scale left in the calculation, say the mass parameter m, so that its actual value becomes irrelevant, since it can be used as energy units and is eventually determined by a comparison with the phenomenology. For instance, sharing the same units of the lattice data, a value m = 0.656 GeV was established in previous works [21,24].
At a finite temperature T = 0, there is a third energy scale and the optimal parameters m, µ become two independent functions of temperature, m(T ), µ(T ), since their optimal ratio is expected to depend on T . In principle, one could proceed as for T = 0 and fix the optimal ratio by monitoring the gauge-invariance of the poles. However, that would at least require a knowledge of the thermal propagators in a generic covariant gauge, while the present formalism has been developed only in the Landau gauge. Moreover, no lattice data are available for a comparison in a generic gauge and finite T . This is not a theoretical limitation by itself, but leads to a weakening of the control of the accuracy.
That of the gauge-invariance of the poles actually is an additional problem one encounters when extending the theory to finite T [28][29][30][31]. Even though the poles of the propagator are constrained to be non-perturbatively gauge-independent by e.g. the Nielsen identities [32], in the thermal formalism different powers of the coupling constant coexist at the same loop order when hard-thermal-loop effects are taken into account, so that consistent resummation schemes are needed in order to obtain truly gauge-invariant results for the poles' position. To first order in the coupling, this can be shown to only affect the imaginary part of the dispersion relations, i.e. the gluon's damping rate. In this work no attempt has been made to implement such resummation schemes or to keep under control the accuracy of the approximation with respect to the issue of gauge-invariance. Whereas at low, non-zero temperatures the screening provided by the gluon's mass may somewhat suppress the effects of the required resummed terms, at higher temperatures the latter are expected to become non-negligible, causing our predictions for the gluon damping rate to become less and less reliable as the temperature is increased.
In the Landau gauge, we explored two complementary strategies and checked that the qualitative description which emerges is robust enough and does not depend on the optimization choice. The first, simpler, strategy consists in using the same m and µ parameters that work at T = 0. That choice was already made in Ref. [26] (albeit with different values for the parameters) and makes sense at low temperature where we expect that m(T ) ≈ m(0) and µ(T ) ≈ µ(0). With this choice, we find the correct qualitative behavior without any adjustment of parameters. In particular, the longitudinal propagator shows a non-monotonic behavior with a crossover at T /m(0) ≈ 0.15. However, the agreement with the lattice data is not quantitative, and the predicted transition temperature is too small (T ≈ 100 MeV), thus indicating that we are already outside the safe low-temperature range. Nonetheless, the disagreement can be absorbed in part by a temperature-dependent optimization of the expansion.
Thus, as a second strategy, we relax the constraints of m and µ being equal to their T = 0 values and regard m(T ) and µ(T ) as independent unknown functions. Reversing the argument that led to their optimization at T = 0, we tune the unknown functions in the Euclidean space by looking for the best agreement with the lattice data. Then, assuming that the higher-order terms are smaller when the agreement is better, the optimized propagators are continued to Minkowski space where the pole location gives information on the dispersion relations of the quasi-gluons at finite temperature. We anticipate that, from a strictly quantitative point of view, the agreement with the lattice is not comparable with the excellent result which was reached at T = 0. Moreover, while the transverse propagator is generally well described, the longitudinal projection becomes very poor deep in the IR for moderately high temperatures. Since most of the deviation occurs below 500-700 MeV, we expect that the predictions for the pole position at high momenta might not be affected too much. We stress that there are no data available in the Minkowski space for a comparison, thus evidencing the power of the method for exploring the analytic properties of the propagators.
Irrespective of the optimization criterion, we confirm the finding of Ref. [26] and the quasi-gluon scenario which was described by Stingl [33], with a gluon which has a very short finite lifetime and can only exist as a shortlived intermediate state at the origin of a gluon-jet event. This paper is organized as follows. In Sec. II we review the set-up and main features of the screened massive expansion and its extension to finite temperatures. In Sec. III we present our results for the Landau gauge gluon propagator at T = 0 and vanishing Matsubara frequency, ω = 0. In Sec. IV we derive the dispersion relations for the quasi-gluons at finite temperatures. In Sec. V we discuss our results and present our conclusions. In the Appendix we explicitly compute the gluon polarization and ghost self-energy at finite temperatures using the screened massive expansion.

II. THE SCREENED EXPANSION AND ITS EXTENSION TO FINITE TEMPERATURE
In a linear covariant ξ-gauge, the gauge-fixed BRST invariant Lagrangian of pure Yang-Mills SU(N) theory is where and L F P is the ghost term arising from the Faddeev-Popov (FP) determinant. The tensor operator is defined asF where the gauge field operators satisfy the SU (N ) alge-braÂ In the standard PT formalism, the total action is split as S tot = S 0 +S I , where the quadratic part can be written as while the interaction contains three vertices In Eq. (5), the standard free-particle propagators for gluons and ghosts, ∆ 0 and G 0 respectively, are defined by their Fourier transforms where the transverse and longitudinal projectors are used Later, we will take the limit ξ → 0 and use the Landau gauge which is a Renormalization Group (RG) fixed point and is the most studied gauge on the lattice. In the above equations, the fields and the coupling must be regarded as renormalized objects and the inclusion of the usual set of counterterms is understood in the total Lagrangian. The massive screened version of PT was developed in Refs. [18][19][20]. At T = 0 and in a generic covariant gauge, the method is very accurate and predictive if the expansion is optimized by the constraints of BRST symmetry [21,24,25]. The expansion arises by a mere change of the expansion point of ordinary PT. Following Refs. [19,21], the new massive expansion is recovered by just adding a transverse mass term to the quadratic part of the action and subtracting it again from the interaction, leaving the total action unchanged. In more detail, we add and subtract the action term where the vertex function δΓ is a shift of the inverse propagator, Figure 1. Two-point graphs with no more than three vertices and no more than one loop. The cross is the transverse mass counterterm of Eq. (13) and is regarded as a two-point vertex.
In the Appendix, a detailed description of the calculation at finite T is given for all the polarization graphs in the figure. and ∆ m µν is a new massive free-particle propagator, Adding that term is equivalent to substituting the new massive propagator ∆ m µν for the old massless one ∆ 0 µν in the quadratic part. Thus, the new expansion point is a massive free-particle propagator for the gluon, which is much closer to the exact propagator in the IR. The mass-shift parameter m is irrelevant in the UV, but acts as a natural cutoff which screens the theory in the IR. Of course, in order to leave the total action unaffected by the change, the same term is subtracted from the interaction, providing a new interaction vertex −δΓ, a twopoint vertex which can be regarded as a new counterterm. Dropping all color indices in the diagonal matrices and inserting Eq. (8) and (12) in Eq. (11), the vertex is just the transverse mass shift of the quadratic part, and must be added to the standard set of vertices arising from Eq. (7). The new vertex is now part of the interaction, even if it does not depend on the coupling. Thus, the expansion has the nature of a δ-expansion, since different powers of the coupling coexist at each order in powers of the total interaction. The proper gluon polarization and ghost self energy can be evaluated, order by order, by the modified PT. In all Feynman graphs, any internal gluon line is a massive free-particle propagator ∆ m µν and the new insertions of the (transverse) two-point vertex δΓ µν are denoted by a cross, as shown in Fig. 1. For further details we refer to Refs. [18,19,21].
Since the total gauge-fixed FP Lagrangian is not modified and because of BRST invariance, the longitudinal polarization is known exactly and is zero. At T = 0, the exact polarization and the dressed gluon propagator are defined by a single function, so that, in the Landau gauge, the exact gluon propagator is transverse, and defined by the scalar function ∆(p). This feature is lost at any finite temperature T > 0, since Lorentzinvariance is broken, and two scalar functions are required instead. In that perspective, it is convenient to maintain the Lorentz structure explicit and to switch to the Euclidean formalism. Then, denoting with p 2 the Euclidean squared momentum, the exact (dressed) gluon and ghost propagators can be written as where t µν and ℓ µν are the Euclidean projectors of Eq. (A5). The proper gluon polarization Π µν and the ghost self-energy Σ are the sum of all one-particleirreducible (1PI) graphs in the screened expansion, including all counterterms. In Fig. 1, the two-point 1PI graphs are shown up to one-loop and third order in the delta-expansion. In the exact self energies, we can single out the tree-level terms and write where the first term m 2 t µν (p) is the tree graph (1a) in Fig. 1 and arises from the insertion of the new two-point vertex −δΓ µν of Eq. (13). We observe that this first tree term cancels the mass shift of the gluon propagator in Eq. (16). Indeed, the physical mass of the gluon arises from the loops and is not merely given by the mass-shift parameter m 2 . The other tree-level terms, −p 2 t µν δZ A , p 2 δZ c , are not shown in Fig. 1 and are the usual field-strength renormalization counterterms. Their UV diverging parts are not affected by the mass parameter and are the same of standard PT [18,19]. The proper functions, Π loop µν , Σ loop , are given by the sum of all 1PI graphs containing loops. The finite parts of δZ A , δZ c are arbitrary and depend on the scheme and on the renormalization scale µ [24,25]. The diverging parts of δZ A , δZ c cancel the UV divergences of the functions Π loop µν /p 2 and Σ loop /p 2 which become finite dimensionless functions of the variable p µ /m. They are defined up to a constant which depends on the dimensionless renormalization scale parameter t = µ 2 /m 2 . Thus, at T = 0, there are two energy scales in the calculation, m and µ. For instance, in a momentum subtraction scheme (MOM) and in the Landau gauge, the one-loop dressed propagators can be written as having made explicit the dependence on N and g 2 as factors in the one-loop functions Π (1) , Σ (1) , according to the notation of Appendix A, where all details of the calculation are reported. In Eq. (18), an explicit choice has been made for the finite parts of the renormalization constants δZ A , δZ c . Of course, that choice depends on the scheme and on the renormalization scale µ. A more general way to get rid of all the scheme-dependent parameters, including the renormalized coupling g 2 , was discussed in previous papers on the screened expansion [18,19,21,24], where two dimensionless one-loop functions were defined (see Appendix B.1 for their explicit expressions), so that the one-loop propagators in Eq. (18) can be recast as functions of the dimensionless variable s = p 2 /m 2 , where z π and z σ are irrelevant normalization constants while all the scheme-dependent parameters are embedded in the two constants π 0 and σ 0 . With some abuse of language, we will refer to them as renormalization constants. Eq. (20) is quite general since it does not require any specific renormalization scheme to be defined. Of course, our ignorance about those constants reflects a well known weakness of the one-loop approximation which depends on the details of the renormalization scheme and on the actual value of the renormalization scale µ. In this sense, we still have two scales, m and µ, and the arbitrary choice of their ratio t = µ 2 /m 2 somehow determines the actual value of the renormalization constants π 0 and σ 0 . A nice feature of the one-loop result is its apparent tangency to the exact result which is approached for special values of the renormalization constants. Those values are equivalent to a choice of the best renormalization scale µ, where the approximation is more effective. It is just an example of the optimized perturbation theory by variation of the renormalization scheme [34,35]. There might be a special scale µ where the expansion converges more quickly and the higher order terms are minimal. Thus, from first principles, we could determine the optimal constants by monitoring some identities which must be satisfied by the exact propagators. For instance, in Ref. [21], the Nielsen identities [36,37] were used, which are a direct consequence of BRST symmetry. From the identities, one can prove the gauge-parameter-independence of the poles and residues of the exact gluon propagator [21]. Then, we might expect that the renormalization constants are optimal when the poles have a minimal sensitivity to the gauge parameter. It is remarkable that the optimized one-loop propagators turn out to be in excellent agreement with the lattice data in the IR. Notably, while the comparison with the data requires an analytic continuation to the Euclidean space, the poles are found in the complex plane. Thus, the one-loop propagators in Eq. (20) are not just one of the many interpolation formula for the data, but they provide a very accurate analytic function in the whole complex plane. The existence of complex poles is one of the most important predictions of the screened expansion. While a thermal mass and a finite damping rate are expected by PT at high temperature, the existence of finite intrinsic values at T = 0 can be regarded as a proof of confinement as first discussed by Stingl [33]. The quasi-gluon has a finite lifetime and can only exist as a short-lived intermediate state. However, at finite temperature, the quasi-gluons play an important role for determining the thermal properties of the hot plasma. Thus, a finite temperature extension of the screened expansion is required for a full study of the dispersion relations which emerge from the pole location.
At a finite temperature T > 0, Eqs. (16), (17) are still valid, but the one-loop graphs in Fig. 1 acquire a finite thermal part which must be added to the vacuum (diverging) contribution at T = 0. The thermal parts are finite and no further renormalization is required. We only have to add the thermal parts to the self-energies in Eq. (17).
We write the Euclidean four-vector as p µ = (p, ω) where ω = p 4 = −ip 0 , while the Lorentz four-vector was (p 0 , p). In the finite-temperature formalism, ω = ω n = 2πnT and the Euclidean integral is replaced by a sum over n and by a three-dimensional integration, Since Lorentz invariance is obviously broken, we introduce a transverse projector P T µν , orthogonal to the fourth Euclidean direction, and its longitudinal complement P L µν , as defined in Eq. (A4), so that the gluon polarization and propagator in Eqs. (16), (17) can be written in the Landau gauge, ξ = 0, as where the projected one-loop dressed functions are and Π (1) L,T are the one-loop projected polarizations, evaluated by projection of the one-loop graphs in Fig 1, omitting the tree graphs. As discussed in Appendix B, each graph contributing to Π (1) L,T can be split as where the vacuum part Π Thus, we can generalize Eqs. (19), (20) and define dimensionless functions so that the projections of the one-loop propagator can be recast as In this form Eq. (26) is quite general since it does not require any specific renormalization scheme to be defined. All the scheme-dependent parameters are embedded in the renormalization constant π 0 . It is not obvious that the same scale µ and constant π 0 which were optimal at T = 0 are still optimal at finite T . Indeed, they might depend on T and even take a different value for the different projections. Moreover, the mass parameter m, which was the only energy scale left after optimization at T = 0, might take a value m(T ) which depends on T . Thus we have three energy scales: the optimal µ(T ), the mass parameter m(T ) and T itself. In other words, according to Eq. (26), at any T and in units of m(0) we have two free parameters, the ratio m(T )/m(0) and the optimal renormalization constant π 0 (T ). Having the role of variational parameters, to be optimized, their best values might be different for the two projections.
While at T = 0 the optimal constant π 0 was determined from first principles [21], by requiring a minimal sensitivity of the poles to any change of the gauge parameter, here we have the less ambitious aim of exploring if a set of optimal parameters does exist such that the screened expansion is able to describe the lattice data with reasonable accuracy. Thus, we work in the Landau gauge and, for each value of T > 0, we fix the parameters by a fit of the available lattice data in the Euclidean space.
At low temperature, as we said, we also explored the alternative of maintaining the parameters fixed at their optimal value for T = 0, in order to give a general description at finite T from first principles, without any input from the lattice and from the known phenomenology. Of course, this approach can only be reliable if T is very low and the thermal effects are small. However, even extrapolating at higher temperatures, the qualitative predictions turn out to be in agreement with the data. Thus, the screened expansion is able to capture the main features of gluon thermodynamics at finite temperature. This is a very important aspect, since our final aim will be to extract some dynamical properties of the quasi-gluons, like the dispersion relations, which cannot be measured on the lattice. Moreover, even qualitative properties, like the existence of complex poles, are of central interest for understanding the behavior of the gluon plasma at high temperature and its phase transition.
In order to fulfill that program, once optimized by one of the two alternatives discussed above, the gluon propagator must be continued to the complex plane. This is a straightforward step if the one-loop graphs are expressed as analytic functions of the Euclidean momentum. A very detailed but tedious analytical evaluation of the integrals is reported in the Appendix. Most of the integrals were encountered in a study of the Curci-Ferrari model [38]. We basically use the same method for decomposing the integrals. However, in the screened expansion there are also some different graphs, namely the crossed graphs in Fig. 1, with one insertion of the mass counterterm.
Their explicit expressions are obtained by a derivative in the Appendix.
Unfortunately, at finite T , not all the multidimensional integrals can be evaluated analytically and an external one-dimensional numerical integration cannot be avoided for almost all the one-loop graphs. Thus, as shown in the Appendix, all the graphs can be written as analytic functions which are defined by integral representations. The remaining integration can be carried out numerically for any complex value of the external momentum, provided that no singularity is encountered along the integration path. Actually, in general, the analytic continuation of integral functions is not trivial. As discussed in Ref. [39], we must check that the external integration on the real axis does not cross any singular point of the logarithmic functions. Otherwise, a modified path must be chosen before the analytic continuation can be undertaken. As shown in Ref. [26], by inspection of the explicit expressions, the existence of singular points on the integration path can be ruled out in the present case. For instance, denoting with Ω = p 0 and p µ = (Ω, p) the external momentum in Minkowski space, the analytic continuation of the thermal integral I αβ (y, −iΩ) is defined by the integral representation of Eq. (B30), where y is the external three vector modulus, y = |p|. We can continue the external energy Ω to the complex plane if there are no singular points on the positive real axis of the integration variable. However, some branch cuts might be present, originating at the singular branch point of the logarithmic function in Eq. (B29) which reads where the complex variable z α is defined as z α = iΩ ± i q 2 + α 2 and ǫ 2 y±q,β = (y ± q) 2 + β 2 . Here α and β are masses equal to 0 or m and q is the integration variable. Assuming the existence of a branch point at q = q 0 on the real axis, the latter must satisfy where the ± signs are independent of each other. Taking a complex energy Ω = Re Ω + i Im Ω with Im Ω > 0, the imaginary part of Eq. (28) gives and substituting back in the real part we obtain which is never satisfied unless Im Ω = β = 0. Thus, if Ω is not real, the branch point q 0 cannot be real and the integral over q, on the real axis, defines an analytic function of Ω. The same argument holds for the other thermal integrals in Appendix B. Thus, we can safely continue the numerical integrals from the Euclidean space (Re Ω = 0, Im Ω > 0) to the whole upper half-plane. Moreover, in the large wavelength limit y → 0, there are no branch points at all because the logarithmic function can be written as L β (z α ; y, q) ≈ log [1 + O(y)] and the argument of the log does not vanish if y is small enough.
Having ruled out the existence of singularities along the integration path, the poles of the gluon propagator and the dispersion relations can be easily extracted numerically in the complex plane by the integral representation of the thermal integrals which are derived in Appendix B.

III. THE GLUON PROPAGATOR AT FINITE T
The longitudinal and transverse projections of the polarization graphs entering in Eq. (26) are decomposed as the sum of more basic Euclidean integrals in Appendix A, for all the one-loop graphs of Fig. 1. The explicit thermal parts of those integrals are presented in Appendix B by integral representations. For any given value of the external three-momentum y = p 2 and Euclidean frequency ω = p 4 = 2πnT , the one-dimensional integrals are evaluated numerically by a simple integration on the real axis and the result is inserted in Eq. (26). We will first explore the projected propagators for π 0 and m fixed at their zero-temperature values which were determined from first principles in Ref. [21]. Then, we will show how their values can be optimized by a comparison with the available lattice data.

A. Expansion optimized at T = 0
In the low-temperature limit, we assume that the optimal renormalization constant π 0 (T ) and mass parameter m(T ) can be replaced by their zero-temperature values π 0 = −0.876 and m(0) = m 0 = 656 MeV, as determined in Ref. [21] by requiring a minimal sensitivity of the pole structure to the gauge parameter. Strictly speaking, in the Landau gauge, that condition fixes π 0 , while m 0 is the only energy scale left and is fixed in order to match the energy units of the lattice data.
Let us first explore the behavior of the gluon propagators as a function of T in the limit ω → 0, where p 2 = p 2 , which is the most studied case on the lattice [11,12]. The longitudinal and transverse propagators are shown in units of m 0 in Fig. 2 and Fig. 3, respectively. The former were multiplicatively renormalized by requiring that with µ 0 /m 0 = 6.098 (corresponding to µ 0 = 4 GeV for m 0 = 656 MeV). We observe that, because of the chosen optimization, in the limit T → 0 the longitudinal and transverse propagators coincide and reproduce the lattice data extremely well [18,19,21,22,24,25], so that the low-temperature limit can be regarded as exact. For reference, in Tab. I we report the physical equivalent of the adimensional temperatures T /m 0 used for the plots. We observe a crossover, in Fig. 2, with the longitudinal propagator which increases in the IR for increasing T below T c ≈ 0.15 · m 0 , but sharply decreases above T c . This non-monotonic behavior is a well known feature which has been reported by several lattice calculations [11,12]. The transverse propagator in Fig. 3, on the other hand,   has a monotonic behavior, decreasing for increasing T , again in qualitative agreement with the known predictions of the lattice. Actually, we cannot expect a quantitative agreement at T ≈ T c or larger values, because we are extrapolating the optimization condition which was valid at T = 0. Thus, the correct qualitative behavior of the propagators at high temperature is an encouraging result. A crude estimate of T c is found by using the zerotemperature value m 0 = 656 MeV for restoring the energy units, yielding at the crossover T c ≈ 100 MeV. This value is quite smaller than the known transition temperature T c ≈ 270 MeV which is measured on the lattice [10][11][12]. The difference might well be the consequence of a sub-optimal choice of the renormalization constant, but it could also arise from a change of the mass parameter with temperature or from the more general failure of PT at high temperature. Thus, it becomes relevant to explore whether a more quantitative agreement might be obtained by a tuning of the free parameters.
B. Optimization by a fit of data at finite T As the temperature increases, our previous assumption, m(T ) = m(0), π 0 (T ) = π 0 (0), becomes less valid.  The lattice data were taken from Ref. [11].
In what follows, we turn to fixing the optimal value of the parameters at T = 0 by a fit of the lattice data of Ref. [11]. Since at non-zero temperatures the projections ∆ L (p, T ) and ∆ T (p, T ) have different behaviors with respect to a change in T , we may expect that the optimal values of the parameters will differ depending on which of the two components of the lattice propagator is used for the fit. This is indeed what we found. Of course, since in the subtracted Lagrangian of the present formalism the gluon mass parameter m 2 (T ) is multiplied by the full four-dimensional transverse projector t µν (p), choosing different mass parameters/scales for the two components of the propagators is not allowed from first principles. This issue will be addressed at the end of this section.
In Figs. 4 and 5 we show, respectively, the longitudinal and transverse components of the gluon propagator at ω = 0 (multiplicatively renormalized at µ 0 = 4 GeV), as functions of the three-dimensional momentum |p| = p 2 , with m(T ) and π 0 (T ) as reported in Tab. II. Such values where obtained by a separate fit of the two components to the lattice data of Ref. [11]; the mass parameters  As we can see, once the parameters are tuned to fit the data, the screened expansion is able to reproduce the lattice propagators quite accurately down to momenta of approximately 0.5 GeV. Moreover, the longitudinal propagator still shows the characteristic non-monotonic behavior with respect to a change in the temperature, increasing at fixed momentum below T = T c ≈ 270 MeV and decreasing above T = T c . Below |p| ≈ 0.5 GeV, the transverse propagator is still in good agreement with the data, while the longitudinal one shows significant deviations, especially at high temperatures. From a numerical point of view, why this is so is exemplified in Fig. 6, where we display the longitudinal propagator for T = 458 MeV and different values of the mass parameter 1 . When tuning the mass parameter m(T ), there is a tension between the low-and intermediate-momentum behavior of the propagator: at lower values of m, the propagator is enhanced (resp. suppressed) below (resp. above) |p| ≈ 1 GeV, so that achieving a good match at low momenta results in a loss of accuracy at intermediate momenta. This behavior is actually shared by both the components of the propagator and at every T = 0, albeit being less significant for the transverse component and at low temperatures. In particular, already at T = 458 MeV the optimal longitudinal values of the mass parameter and of the renormalization constant strongly depend on the choice of a lower cutoff momentum for the fit to the lattice data; for this reason, we do not report them.
As anticipated earlier, the optimal mass parameters (and renormalization constants) needed to reproduce the lattice data differ for the two components of the propagator. In Fig. 7 At T = 260 MeV ≈ T c , the optimal value of m(T ) is nearly equal for both the projections, namely m(T ) = 425 − 450 MeV. As for the renormalization constant, except for the point at T = 290 MeV ≈ T c , the optimal π 0 (T ) increases with the temperature when fitted from the transverse propagator. When optimized by the longitudinal propagator, on the other hand, it shows a nonmonotonic behavior, decreasing below T c and increasing again above T c .
The large differences in the optimal values of m(T ) and π 0 (T ) obtained for the two projections make it clear that, in the present formalism, it is not possible to quantitatively recover both the longitudinal and the transverse component of the gluon propagator by a unique choice of parameters. Thus at T = 0 the screened expansion appears to be suboptimal as a "variational" ansatz. At least in part, this could be expected on the basis of what is known about the high-temperature, low-momentum behavior of the Yang-Mills propagators: at large temperatures and low momenta, the gluons' thermal mass is best described by a momentum-and direction-dependent Hard Thermal Loop (HTL) term in the Lagrangian, given by [41] where m 2 el (T ) = g 2 N T 2 /3,ŷ is a light-like four-vector and the integration is over the directions ofŷ. To first order in the coupling, ∆L HTL generates two different thermal masses for the three-dimensional projections ∆ L (p, T ) and ∆ T (p, T ) of the gluon propagator. By not taking into account this difference, the screened expansion lends itself to a breakdown at large temperatures, which can be partially avoided if the mass parameter and renormalization constant are tuned to separately fit the two projections.
The simplest way of solving this issue in the context of the screened expansion, i.e. without resorting to a HTL resummation, would be to change the expansion point of perturbation theory in such a way that the two three-dimensional projections of the zero-order gluon propagator, ∆ T m and ∆ L m , have different masses ab initio. This can be achieved by redefining the kernel δΓ µν (p; T ) = m 2 (T ) t µν (p) of the shift of the action δS as where m T (T ) and m L (T ) are independent massparameter functions for the two projections. With such a prescription, in a general covariant gauge the zero-order Euclidean gluon propagator ∆ µν m (p; T ) would read ∆ m (p; T ) µν → ∆ T m (p; T ) P T µν (p) + ∆ L m (p; T ) P L µν (p)+ where are the sought-after zero-order propagators. Setting-up the perturbation theory with independent mass functions for the two projections would give us the freedom to optimize the former separately from first principles, according to the behavior of the respective dressed propagators. Implementing the shift in Eq. (35), however, is a non-trivial task: having different longitudinal and transverse masses running in the loops breaks the Lorentzinvariance even of the simplest vacuum integrals and, more generally, requires a complete recalculation of the gluon polarization.

IV. DISPERSION RELATIONS AT FINITE T
Being in possession of analytical expressions (modulo a one-dimensional integration at finite T ) for the Euclidean gluon propagator allows us to analytically continue the latter to the whole complex plane so as to study its singularities. As is well known, the location of the poles of the propagator gives us information on the dispersion relations of the gluonic quasi-particles: the energy ε T,L (p, T ) and damping rate γ T,L (p, T ) of the quasi-particles, as functions of the three-dimensional momentum p and of the temperature T , are obtained by solving the equation where ω = ε − iγ (modulo a factor of i) extends the real and discrete Matsubara frequencies ω n = 2πnT to the complex plane and the subscripts T, L refer to the components of the propagator. At non-zero temperatures and momenta, the poles of the two components are expected to be found at different locations, yielding two separate branches of the dispersion relations. The limit T → 0 of the dispersion relations was already studied in the framework of the screened massive expansion in Refs. [20][21][22]. In [21] we found that the zerotemperature gluon propagator (whose longitudinal and transverse three-dimensional components are constrained to be equal by Lorentz simmetry) has two complexconjugate poles at −p 2 = m 2 pole , (m 2 pole ) * , where, setting m 0 = 656 MeV by sharing the same units of the lattice, with m 2 pole = m 2 R + i m 2 I . In terms of ε vac (p) = lim T →0 ε T,L (p, T ) and γ vac (p) = lim T →0 γ T,L (p, T )and singling out one of the poles -, this translates into the dispersion relations Clearly, At the other end of the spectrum, as |p| → ∞, the gluon's vacuum dispersion relations reduce to those of a massless particle, ε vac (p) → |p|, γ vac (p) → 0. Under the assumption that the optimal masses m(T ) and renormalization constants π 0 (T ) reported in the previous section only depend on the temperature, and not on the Matsubara frequency ω n , the finite-T dispersion relations of the gluon quasi-particles can be easily extracted from the screened expansion's gluon propagator, making use of said parameters (cf. Tab. II). We remark that, since at low momenta the longitudinal projection was not found to be in good agreement with the lattice data for any value of the parameters, the longitudinal dispersion relations are expected to be reliable only at sufficiently high momenta (say above |p| ≈ 0.5 − 0.7 GeV).
In Figs. 8 and 9 we plot the energy ε T,L (p, T ) and the damping rate γ T,L (p, T ) of the transverse and longitudinal gluons at fixed T , as functions of the momentum |p|. As we can see, below the critical temperature T c ≈ 270 MeV both the transverse energy and the transverse damping rate (Fig. 8) are suppressed with respect to their zero-temperature (vacuum) limit, with the effect being more pronounced for ε T than for γ T . Above T c this behavior is reversed; the transverse energy starts to approach again its vacuum limit, while the damping rate grows larger than it. The longitudinal branch (Fig. 9) shows a more significant suppression in both the energy and the damping rate below T c , with γ L becoming quite small at high momenta around the critical temperature. At higher temperatures both ε L and γ L start to approach back their vacuum limit. 2 In the limit p → 0 and for any non-zero ω, the longitudinal and the transverse projection of the gluon propagator are known to collapse to a single temperaturedependent function; as a consequence, the corresponding branches of the dispersion relations share the same zeromomentum limit. The p = 0 poles of the gluon propagator are located at −i(ε 0 (T ) − i γ 0 (T )), where are, respectively, the mass and the (zero-momentum) damping rate of the gluon quasi-particles. With regards to such a constraint, the optimized framework of Sec. IIIB is inconsistent: using different mass parameters for the longitudinal and the transverse projections of the propagator causes the two branches of the dispersion relations to have unequal p → 0 limits. All the same, as previously discussed, the low-momentum limit of the longitudinal gluon propagator was found to be quantitatively unreliable at temperatures which are not vanishingly small. It follows that the p → 0 limit of the longitudinal dispersion relations cannot be trusted regardless of the inconsistency. Since only the screened expansion's transverse propagator, with the parameters in Tab. II, was found to reproduce the lattice data at low momenta, in what follows we will make use of the transverse dispersion relations to study the behavior of ε 0 (T ) and γ 0 (T ). From first principles, it is understood that a good description of the long-wavelength longitudinal gluon excitations must yield the same results.
In Fig. 10 we display the mass and the zero-momentum damping rate of the gluon quasi-particles as functions of the temperature. Across the critical temperature, both of them show a characteristic behavior, decreasing below T c and increasing again in a linear fashion above T c . The mass decreases from ε 0 (0) = ε vac (0) = 581 MeV to ε 0 (T c ) ≈ 450 MeV, whereas the zero-momentum damping rate slightly decreases from γ 0 (0) = γ vac (0) = 375 MeV to about 350 MeV around T c . The increase in the damping rate actually seems to start somewhat below the critical temperature (see the data point T = 260 MeV in Fig. 10); we could not determine whether this is a physically meaningful behavior or an artifact due to uncertainties in the parameters of Tab. II. The behavior of the gluon mass in Fig. 10 confirms the picture of a confined gluon -whose mass is dynamically generated through the strong interactions themselves like in the T → 0 limit -which becomes deconfined above the critical temperature T c ≈ 270 MeV. In the deconfined phase, the mass of the gluon is thermal in nature and increases linearly with the temperature. The same qualitative behavior was observed in [26], where the gluon mass and zero-momentum damping rate were studied in the screened expansion at finite T using the same scheme of Sec. IIIA, i.e. taking temperature-independent values for both the gluon mass parameter m and the renormal-ization constant π 0 .

V. DISCUSSION
The comparison with the available lattice data showed that the screened expansion gives a correct qualitative description of the gluon propagator at finite T . The agreement improves if the renormalization constants are tuned at each value of the temperature. At high temperatures and deep in the IR, the failure to reproduce the longitudinal projection might arise from the combined effect of several issues like the need of some HTL resummation, a poor optimization and the inadequacy of the single-mass splitting of the action at a finite temperature. Indeed, the lattice data seem to suggest that a two-mass scheme should be introduced from the beginning for extending the screened expansion at a finite temperature. Nonetheless, the qualitative behavior of the propagators seems to be correct and quite robust, irrespective of the optimization scheme. The pole trajectories can be determined in the complex plane, yielding valuable predictions which cannot be extracted from the lattice data in the Euclidean space. We have reported in some detail the dispersion relations of the quasi-gluon for several temperatures across the deconfinement transition.
An important feature which emerges from our study is a crossover at the deconfinement transition. The energy of the quasi-particle is suppressed by temperature in the confined phase. On the other hand, above the critical temperature, the behavior is reversed and the energy increases as a function of temperature. The same effect can be observed for the physical mass, defined as the long wavelength limit ε 0 (T ) of the pole's real part, as shown in Fig. 10. In the confined phase, the mass decreases like an order parameter being suppressed by the temperature. This behavior is consistent with that of a dynamical mass which is related to a condensate, the latter being expected to vanish at the transition temperature. However, at finite temperature the quasi-gluon is also expected to acquire a thermal mass which increases linearly, like any other quasi-particle. The two effects might coexist across the transition, yielding a crossover rather than a sharp transition. In the low-temperature limit the dynamical nature of the mass dominates, while above the deconfinement transition the mass becomes a pure thermal mass. Thus, we argue that in the low-temperature phase the mass suppression might be a signature of the dynamical nature of the gluon mass. On the other hand, as discussed in Ref. [26], the existence of an intrinsic damping rate, which saturates at a finite value at T = 0, is a confirmation of the quasi-gluon scenario laid out by Stingl [33]. The massive gluon also has a very short finite lifetime and is canceled from the asymptotic states [26], suggesting that the gluon quasi-particles of the interacting vacuum can only travel the short distance of about a Fermi and can only exist as intermediate states at the origin of a gluon-jet event.

ACKNOWLEDGMENTS
The authors are in debt to Orlando Oliveira for sharing the lattice data of Ref. [11].
This research was supported in part by "Piano per la Ricerca di Ateneo 2017/2020 -Linea di intervento 2" of the University of Catania.
Appendix A: one-loop graphs

Notation
The Euclidean four-vector p µ is defined as where ω = p 4 = −ip 0 and the the Lorentz four-vector is (p 0 , p).
In the finite temperature formalism, ω = ω n = 2πnT and the Euclidean integral is replaced by a sum over n and by a three-dimensional integration The generic (massive) propagator G m (p) is At finite temperature, it is useful to introduce the following orthogonal projectors beside the Lorentz projectors The trace of the projectors is The dressed Euclidean propagator of the gluon can be written as ∆ ab µν (p) = δ ab ∆ µν (p) where and the gluon polarization is Π ab µν (p) = N g 2 δ ab Π µν (p). Since Π µν (p) is transverse, i.e. p µ Π µν (p) = 0, in the Landau gauge (ξ → 0) the dressed propagator is also transverse. We introduce the projected polarizations so that the total polarization reads and the dressed propagator can be written as where the projected parts are In the Landau gauge, ξ → 0, the propagator is transverse and its components are determined by the projected polarizations Π T (p) and Π L (p). The graphs are evaluated in the Landau gauge, using the (transverse) massive free propagator [G m (p) t µν (p)] in the internal gluon lines. The dressed Euclidean propagator of the ghost can be written as G ab (p) = δ ab G(p), where and the ghost self energy is Σ ab (p) = δ ab Σ(p). In the graphs, the massless free propagator −G 0 (p) is used in the internal ghost lines. All the uncrossed one-loop graphs can be decoupled by the method of Ref. [38] and written in terms of the set of integrals together with their projections I αβ L (p) = P L µν (p) I αβ µν (p), Explicit expressions are reported in Appendix B. By exchanging k µ and p µ −k µ in the integrals, it is easy to show that I αβ (p) = I βα (p), while in general I αβ µν (p) = I βα µν (p). However, since p µ P L,T µν (p) = 0, and the projected integrals turn out to be symmetric, I αβ L,T = I βα L,T . We note that I αβ L p and I αβ T p might depend on p because of the explicit dependence in the projectors. For instance, let us consider any constant integral which does not depend on the external momentum p. Let us denote by I L,0 , I T,0 , the non-zero components that can be written, taking k µ = (k, ω n ), as In fact, the explicit projections I L,p , I T,p can be defined and evaluated as in Eqs. (A14): More generally, for the integral I αβ µν (p), which has an explicit dependence on p, the projections have the following limits: where I αβ L,0 = I αβ 44 (0), I αβ T,0 = I αβ ii (0).
The limits in Eq. (A20) agree with the physical requirement that transverse and longitudinal projections must coincide for any ω in the limit p → 0, while they are different for any finite p in the limit ω → 0. Each crossed graph Π × , containing one insertion of the mass counterterm, can be obtained by the corresponding uncrossed graph Π by a simple derivative Their explicit calculation requires the definition of a new set of integrals ∂I αβ , ∂I αβ L,T , ∂J m , ∂ 2 J m : We note that the second argument (β) is kept fixed in the derivative, so that ∂I αβ = ∂I βα . When α = β the derivative must be taken twice, so that for instance (A25) Not all the integrals are independent. For instance, it can be easily shown that It is useful to introduce the integrals J L,T m which follow by setting f (k) = G m (k) in Eq. (A17), so that Eqs. (A26) can be extended to the projected integrals, and, by Eq. (A21), the projections I αβ L,T p can be expressed in terms of the constant integrals J L,T m .

Graph 1b -(tadpole)
Setting d = 4, Eq. (31) of Ref. [40] reads yielding where the integrals J α , I αβ µν (p) were defined in Eqs. (A13) and their explicit expressions are reported in Appendix B. The projected polarization of graph (1b) is The vacuum contribution can be extracted by evaluating the integrals in the limit T → 0 where I m0 µν (0) → 1 4 δ µν J m so that in agreement with the general result of Ref. [19] for d = 4.

Graph 2b -(gluon loop)
The general explicit expression for the graph (2b) has been reported in Ref. [40], for a generic dimension d and a generic free-particle propagator. In the Landau gauge, the explicit expression for d = 4 can be written as (see also Ref. [38]) where, denoting q = p − k, All integrals can be evaluated by the method of Ref. [38] and written in terms of the integrals in Eq. (A13).
Adding Eq. (A45) and Eq. (A48) in Eq. (A41), the second polarization term in Eq. (A34) is * * * * * * * The third polarization term in Eq. (A34) can be decomposed by observing that so that, changing the integration variable from k to q in the q ν p µ term, Π 3 reads The first integral can be decomposed by using the identity and observing that, by the second of Eqs. (A36) and the second of Eqs. (A43), yielding In the second integral Π 3b we can use the identity where the last term can be dropped because it is antisymmetric in the exchange of k and q and its contribution to the integral is zero. Taking the second of Eqs. (A43) with α = β = m and the same equation with α = β = 0 and k,q interchanged, their product can be written as where the second of Eqs. (A36) has been used for decomposing the products of more than two G functions. Then, the integral can be written Adding Eq. (A 3) and Eq. (A58) in Eq. (A51), the third polarization term in Eq. (A34) is * * * * * * * The last polarization term in Eq. (A34) can be decomposed by observing that Then, recalling that q µ = p µ − k µ , the integral reads Using the identity the two pieces can be added together yielding The product of three G functions can be decomposed by the second of Eqs. (A36) and the two arising terms can be written by the first of Eqs. (A36), with (α, β) = (m, 0) and (α, β) = (m, m), respectively, The integral then reads where Using the first of Eqs. (A43) with α = 0 and β = m, 0, and decoupling the product G m (k)G 0 (k) by the second of Eqs. (A36), the term Π 4a can be written as so that the integral reads Π 4aµν (p) = p µ p ν p 4 2m 4 I 00 (p) + The transverse projections of the graph follow by the projected integrals in Eqs. (A14). We observed that by Eq. (A15) the projected integrals turn out to be symmetric, I αβ L,T = I βα L,T . Thus, the projection of the graph follows by dropping the longitudinal and the antisymmetric terms, and by replacing the integrals by the projected ones according to

Graph 2a -(ghost loop)
In the Landau gauge, setting d = 4 and using a freeparticle propagator, the general expression of the ghost loop (see e.g. Ref. [40]) reads By exchanging k µ and p µ − k µ the integral shows the symmetry Π µν = Π νµ so that, using Eq. (A53), we can replace (A74) The first two terms on the right hand side cancel in the integration yielding The projected ghost loop is just

Ghost self-energy
In this work, the total one-loop ghost self energy is the sum of the standard one-loop graph and the crossed one, which contains the insertion of a mass counterterm, where Σ(p) is the standard one loop integral [19,40] in the Landau gauge, Using the last of Eqs. (A43) with α = m and β = 0, and decoupling the product G m (k)G 0 (k) by the second of Eqs. (A36), we can write Then using the second of Eqs. (A43) with α = β = 0 and dropping the vanishing integrals the second term of Eq. (A79) reads and the (uncrossed) one-loop self energy can be written as as derived in Ref. [38] by the same method.
Finally, the doubly crossed tadpole (1d) in Eq. (A84) can be written as [19,40]  is the Euclidean integral at T = 0, denoted vacuum part I V = I(0), while the thermal part I T h (T ) is Resid.
2 ℜf (k, ik 0 ) e βk0 − 1 Rek0>0 (B3) where the sum is over the residues in the right complex plane of k 0 and the symbol ℜf is defined as We observe that if f (k) is a complex function, then ℜf (k) is not the true real part Ref (k). The thermal part vanishes in the limit T → 0.
Many of the thermal integrals were evaluated in great detail in Ref. [38]. In the next sections we collect the same results and, by the same method, we add the explicit evaluation of all the remaining integrals that are required in the present work.