Non-perturbative finite-temperature Yang-Mills theory

We present non-perturbative correlation functions in Landau-gauge Yang-Mills theory at finite temperature. The results are obtained from the functional renormalisation group within a self-consistent approximation scheme. In particular, we compute the magnetic and electric components of the gluon propagator, and the three- and four-gluon vertices. We also show the ghost propagator and the ghost-gluon vertex at finite temperature. Our results for the propagators are confronted with lattice simulations and our Debye mass is compared to hard thermal loop perturbation theory.

In this work, we focus on thermal correlation functions of the pure gauge sector of QCD. Quantitative control over Yang-Mills theory at zero as well as finite temperature is a pivotal prerequisite for predictive investigations of the QCD phase structure with functional methods. While vacuum Yang-Mills correlation functions have been studied intensively in the past two decades [2,, results for the finite-temperature correlation functions are scarce, see [61][62][63][64][65][66][67][68][69][70][71] for propagator studies. For the vertices, the situation is even less satisfactory and only exploratory studies exist [72,73].
The main goal of this study is to get quantitative access to the finite-temperature 1PI n-point functions of Yang-Mills theory. These correlators contain all the information about the observables. For example, the resulting propagators and vertices can be used to investigate the center-symmetry phase transition in terms of the Polyakov-loop potential, see e.g. [74][75][76][77][78][79][80][81][82][83][84][85]. Furthermore, the Debye mass, which has been studied intensively on the lattice [86][87][88][89] and hard thermal loop perturbation theory [90,91] as well as with other thermal QCD approaches [92][93][94], can be extracted from the gluon propagator. Additionally, the correlators allow for the extraction of spectral functions and the calculation of the shear viscosity [95,96].
To calculate the 1PI n-point functions, we perform a systematic vertex expansion of the effective action with the aim of quantitative precision, controlled by apparent convergence. The zero-temperature baseline for this calculation is provided by [2], a recent FRG study, which incorporates all tensor structures present at the classical level in a self-consistent truncation scheme. Here, we generalise this truncation to finite temperature, which includes the splitting of the correlation functions into electric and magnetic components. In particular, we provide results for the electric and magnetic gluon propagators as well as the electric and magnetic components of the threeand four-gluon vertices. For the propagators, we compare extensively to results obtained in lattice simulations. We use the Debye screening mass to determine a lower bound for the temperature range in which hard thermal loop perturbation theory can be applied straightforwardly. Furthermore, the finite-temperature behaviour of the ghost-induced zero crossing of the three-gluon vertex is investigated. The comprehensive truncation brings along new technical challenges whose solutions are discussed. In summary, this work provides a major step towards investigations of the QCD phase structure from first principles within the functional methods.
This paper is organised as follows. In Sec. II we discuss the finite-temperature vertex expansion, order parameters, and the Debye screening mass. Section III deals with finite-temperature flows of gauge theories. We present our results in Sec. IV and discuss them in Sec. V. Finally, we summarise our findings and give an outlook in Sec. VI. Technical details and numerical checks are provided in the appendices. In particular, we confirm regulator independence in App. A.

II. YANG-MILLS THEORY AT T > 0
We consider Euclidean Yang-Mills theory, whose classical action in general covariant gauges is given by Here, A, c andc denote the gluon, ghost and antighost fields and x = d 4 x . The gauge fixing parameter ξ is taken to zero in Landau gauge, ξ → 0 . The field strength tensor and adjoint covariant derivative are given by where f abc are the structure constants of the Lie algebra. Our notation largely follows earlier works within the fQCD collaboration [23], and we refer to [1][2][3]26] for further details.

A. Finite-temperature vertex expansion
Functional approaches require an approximation of the corresponding generating functional. We use a vertex expansion about the vanishing expectation values of the gluon and ghosts fields, A µ = 0 and c =c = 0 . These field values are solutions of the equations of motion and constitute the vacuum at vanishing temperature. The intricacies at finite temperature are discussed in more detail in the next Sec. II B. In the vertex expansion, the effective action is written as a sum over powers of the fields, Γ[Φ] = n 1 n! {pi} Γ (n) (p 1 , . . . , p n ) Φ(p 1 ) · · · Φ(p n ) , (3) where Φ = (A µ ,c, c) is a superfield and momentum conservation implies i p i = 0 . The expansion coefficients in (3) are the 1PI n-point functions that are in field components given by Γ (n) The correlation functions are expanded in terms of basis tensors T (i) and dressing functions λ At finite temperature, the vacuum O(4)-symmetry is replaced by Z 2 ⊗ O(3) . This reduced symmetry implies a difference between the magnetic and electric components, which correspond to the directions that are transverse and longitudinal with respect to the heat bath. Starting FIG. 1: Constituents of our vertex expansion. We use the classical tensors that are present in the bare action and attach magnetic (blue) and electric (red) projection operators to the gluon legs. Missing combinations, e.g. vertices with one electric leg, vanish if the Matsubara modes are set to zero and are not computed in our truncation. .
from the longitudinal and transverse vacuum projection operators, we decompose four-vectors into where n ∈ Z are the discrete Matsubara modes and ω n = 2π T n the corresponding frequencies. This leads to the magnetic and electric projection operators at finite temperature, A crucial consequence of the breaking of the vacuum O(4)-symmetry by (8) is the splitting of the tensor structures into electric and magnetic components. In particular, the propagators are given by with dimensionless scalar dressing functions 1/Z M A and 1/Z E A for the magnetic and electric components of the gluon propagator. In the case of the vertices, we take only the classical tensor structures into account. Similarly to the gluon propagator, we split their dressings into electric and magnetic components. See Fig. 1 for an illustration of the constituents of our truncation and App. B and C for further details. As a consequence of the restriction to classical tensors only, the tensor bases of the gluonic vertices are not complete, and the projection of the tensor equations onto the scalar dressing functions is not unique. We use vacuum calculations to identify uncertainties that stem from the projections in order to disentangle them from finite-temperature effects, see App. D for details.
Due to the breaking of O(4)-invariance, the scalar dressings are in general functions of the Matsubara modes and spatial momenta, e.g. Z(p) = Z(ω 2 n , p 2 ) for a generic wave function renormalisation Z. However, the thermal contributions to the correlation functions decay rapidly for spatial momenta and frequencies with p 2 > ∼ (2π T ) 2 . Hence, the thermal correlation functions converge quickly towards their O(4)-symmetric vacuum counterparts for these momenta, see e.g. Fig. 6 and 9. Consequently, the correlation functions exhibit an approximate O(4)-symmetry for all higher Matsubara modes, and most of the finite temperature effects are encoded in the zero mode at small spacial momenta p 2 < ∼ (2π T ) 2 . Therefore, the spatial momentum dependence of the Matsubara zero modes can be used to obtain a very good approximation of the full frequency and momentum dependence via or in short Z(p) = Z(0, p 2 ) . In this work, we compute the zero modes of the propagators and use (10) to close the equations. Within functional methods, this O(4)symmetric approximation has been found to be quantitatively reliable for gluon [62,97] as well as quark propagators [67]. This is confirmed by lattice studies that show a slight deviation of the O(4) invariance only for the first Matsubara mode at temperatures just below the critical temperature [98,99]. This pattern carries over to the scalar dressings of higher order correlation functions λ (n) (p 1 , . . . , p n ) . Analogously to the propagator dressings, we base our computation on the zero modes λ (n) ( p 1 , . . . , p i ) = λ (n) (ω n1 = 0, p 1 , . . . , ω nn = 0, p n ) .
In contradistinction to the propagator dressings, the zero modes of the vertex dressings λ (n) depend on all p i · p j and not only p 2 . However, the spatial momentum dependence of the vertices is well described by a onedimensional symmetric-point approximation in four [2] as well as three dimensions [100,101], see also App. A. This leads to which allows to compute the flows of the zero modes of the vertex dressings in a quantitatively reliable approximation, cf. Fig. 12. However, the flows of the zero modes depend on the full frequency and spatial momentum dependence. Analogously to the propagator dressings, we approximate the full momentum dependence with an O(4)-symmetric generalisation of (12), where the symmetric momentump is then given bȳ In summary, we use two quantitatively reliable approximations for the dressing functions: the approximate O(4)-invariance of all non-vanishing Matsubara modes, which allows to use only information from the lowest Matsubara mode, and the well-tested symmetric point approximation. This truncation generalises the vacuum truncation used in [2], see App. E for an explicit numerical check.

B. Non-trivial vacuum and backgrounds
As discussed in the last Sec. II A, we use a vertex expansion about vanishing field expectation values A µ = 0 and c =c = 0 . This necessitates a thorough discussion of the implications of this choice, in particular for comparisons to lattice results. We argue that such an expansion about vanishing background fields, i.e. Landau gauge, leads to correlation functions that agree with the lattice correlators for temperatures outside a small region around the phase transition. Furthermore, even near the phase transition, sizeable effects are only expected for correlation functions that have electric gluon legs, the electric gluon propagator being their most prominent representative. This becomes most evident by investigating the relation of the physical solution of the equation of motion in non-vanishing gluon background fields and the Polyakov loop, the canonical order parameter of the confinement-deconfinement phase transition. For the convenience of the reader, the first two parts briefly review corresponding relevant advances in non-perturbative functional approaches [74][75][76][77][78][79]102], see [65,[80][81][82][83][84][85] for further applications.

Correlation functions
To facilitate the discussion, we use the background extension of Landau gauge, called Landau-deWitt gauge. Here,Ā µ is a general background and a µ is the quantum fluctuation field. In this formulation, the effective action is gauge invariant under background gauge transformations, which allows for a simpler interpretation of physical backgrounds as well as simpler technical implementations. Besides being a functional of Φ = (a µ ,c, c) , the effective action depends now also on the backgroundĀ . Accordingly, the vertices Γ (n) [Ā, Φ] are correlation functions in the background Here, we have introduced the background current, which satisfies J(Ā, Φ) = δΓ/δΦ . The correlation functions in the absence of external sources, J(Ā, Φ) = 0 , are then given by Γ (n) [Ā, Φ EoM ] , where Φ EoM is a solution of the equation of motion in the chosen backgroundĀ , In general, this conditions yields stationary points of the effective action. In particular, the expansion point (Ā, Φ) = 0 satisfies (17), but does not necessarily single out the physical minimum. In contrast, the physical correlators that correspond to scattering amplitudes are obtained at the physical solution of the equation of motion (17), i.e. the minimum of the effective action (Ā, Φ min [Ā]) . This is also the field value about which the vertex expansion is expected to be most stable and converge most rapidly. Furthermore, only an expansion around the physical solution of the equation of motion allows for a direct comparison with correlation functions from lattice simulations, since the latter are measured on the physical ground state. In general, any other expansion point requires information about higher correlation functions in order to evaluate Γ (n) [Ā, Φ min ] . In particular, in a vertex expansion with expansion point (Ā, Φ) = 0 , the inverse propagator is given by where we suppressed external momentum arguments. Therefore, we expect deviations between the correlation functions Γ (n) [0, 0] , computed in this work, and those from lattice simulations. However, these differences are sizeable only if the momentum scales of the solution (Ā, Φ min ) = 0 are of the same order as T c , the characteristic scale of the finite temperature Yang-Mills system. Only in this case, the higher correlation functions would lead to noticeable contributions in (18). We can utilise the background field to achieve a technical simplification. Since it is arbitrary, we can choosē For this particular choice, the background carries all the non-trivial information about the ground state, whereas the (classical) fluctuation field vanishes on the equation of motion. The physical correlators are then given by Γ (n) [ Ā , 0] . In particular, the inverse propagator (18) for the gluon is then given by  [79]. Both observables are order parameters for the confinement-deconfinement phase transition. Moreover, Semi-perturbative studies in the Curci-Ferrari model for Yang-Mills theory confirm that the background has large effects on the electric propagator at temperatures close to the phase transition [65]. Furthermore, the calculation of quantitatively correct values for the chiral phase transition temperature as well as its observed coincidence with the confinement-deconfinement crossover temperature require to take into account such a non-trivial minimum [76]. Finally, such a consistent treatment was also required for the description of the Roberge-Weiss transition [76] as well as the study of criticality in SU (2) Yang-Mills theory [75,102].

Order parameters
A further advantage of the background Ā is its relation to the Polyakov loop [74,76,78,79], the standard order parameter of the confinement-deconfinement phase transition. The traced Polyakov loop is expressed as a correlator of the temporal gauge field with  expectation value of the algebra element takes the particularly simple formφ for a given backgroundĀ 0 . Due to background gauge invariance, the eigenvalues ofφ can therefore be calculated from the eigenvalues ofĀ 0 in any background gauge. In particular, the effective potential, V [Ā 0 ], of Landau-deWitt gauge carries thus the full information about the eigenvalues of the expectation value of ϕ .
In conclusion, the effective potential V [Ā 0 ] is an order parameter potential for center symmetry. The gauge invariant observables, L[A 0 ] and L[ Ā 0 ], or equivalently also Ā 0 , serve as order parameters for the confinementdeconfinement phase transition, see Fig. 2. Therefore, the vanishing of L[A 0 ] in the confined phase relates to a non-vanishing value for Ā 0 . This has recently been demonstrated explicitly in a self-consistent vertex expansion scheme, which has been used for the first computation of L[A 0 ] within functional methods [79]. Finally, the electric propagator A 0 (p)A 0 (−p) is closely related to the propagator of an order parameter field, and as such should show critical properties, see [78]. Hence, we expect the electric correlators to be affected most by the background field.

Comparison to lattice simulations
The previous discussion of the non-trivialĀ 0 background and its relation to the order parameter of the confinement-deconfinement phase transition allows us to derive a theoretical estimate of the temperature range, in which our present results potentially deviate from the respective lattice results due to the different background configurations. The first important piece of information is given by the fact that the order parameter L[ Ā 0 ] approaches unity rapidly for temperatures above the phase transition temperature, see Fig. 2. This in stark contrast to the Polyakov loop L[A 0 ] , which is usually calculated in lattice simulations. The latter reaches its asymptotic value only for T T c , which can be understood from fluctuation effects [79]. The fact that L[ Ā 0 ] quickly approaches unity above the transition temperature can be formulated as the more precise statement, As a consequence, we can expect quantitative effects due to the non-trivial background only at temperatures T < ∼ 1.3 T c . The most immediate effect of this nontrivial background is a shift in the Matsubara frequencies ω n → ω n ± 2π T ν i , where ν i are the eigenvalues ofφ , or equivalently of βg Ā 0 /(2π) . Rotating the constant field into the Cartan sub-algebra, these are given by in SU (2) and SU (3) , see e.g. [79]. However, for T < ∼ 0.5 T c the effect of the shifts of the Matsubara frequencies is suppressed by the zero temperature gapping m gap of the gluon propagator 2πν i T /m gap 1 . Therefore, we expect sizeable effects due to the non-trivial background only in the regime and in particular in the electric gluon propagator.

C. Debye screening mass
Gluons are screened at high temperatures by the standard thermal Debye mass. However, also in the con-fined phase, they posses a finite screening mass. Our non-perturbative results allow to compute a screening mass also below the critical temperature. We extract it from the zero mode of the electric gluon propagator, , whose computation is detailed below in Sec. III. To this end, we Fourier transform the propagator, At high temperatures, the screening mass can then be extracted from the exponential decay at large distances, The screening mass m s obtained with (27) is shown in Fig. 3. The large distance behaviour of G E T (x) and the fits by (27) are provided in App. I. The left panel shows that the screening mass is finite across the phase transition and possesses a minimum at some finite temperature. Perturbatively, the Debye mass is given to leading order by A prescription for taking higher-order effects into account has been proposed by [92] In order to compare our screening mass to the expressions (28) and (29), we have to determine g T and c D . We use and fit c D at large temperatures to our result since it is not computable within perturbation theory. The running coupling of the (electric) three-gluon vertex α E A 3 is introduced below in (37).
As shown in Fig. 3b, the Arnold-Yaffe Debye mass agrees almost perfectly with our non-perturbative result down to T ≈ 0.6 GeV. In contrast, the leading-order Debye mass deviates instantly from our result. By default, we set c p = 1 in (30) because this is the scale that is

III. METHOD
In this section, we discuss the flow equations, the implications of the regulator term at non-vanishing RG scale, and we provide details on the numerical implementation.

A. FRG flows
The functional renormalisation group in Wetterich's formulation [104] allows to integrate momentum-shell contributions to the effective action in the Wilsonian spirit. To this end, a scale-and momentum-dependent mass term is added to the classical action, The regulator term suppresses quantum as well as thermal fluctuations at momenta below the RG scale k . Taking the derivative with respect to the scale k leads to [69,103].
where t = ln(k/Λ) denotes the RG time and Using the Matsubara formalism, the momentum integral in (32) is given by where q 0 = 2π T n ≡ ω n . A graphical representation of the Wetterich equation for the effective action is shown in Fig. 4. The truncated flow equations for the correlation functions that are obtained by taking functional derivatives of (32) are displayed in Fig. 5. Instead of the flat regulator [105] used in the vacuum computation [2], we use an exponential regulator. As demonstrated in App. A, the results for the correlation functions do not depend on this choice within our error bars. However, analytic regulators such as the exponential regulator are better suited for numerical calculations of thermodynamic quantities since they carry the thermal exponential decay with the cutoff scale ∼ e −ck/T in the flow [62,97], see [106] for a detailed study.
To reduce the numerical effort of the finite-temperature calculation, we exploit the degeneracy of the dressings for k 2π T . We integrate the finite-temperature flow starting from the non-trivial zero-temperature effective action at with λ = 4 and Λ min T = 1 GeV, see App. F for details.

B. Renormalisation and mSTIs
In the presence of a regulator, the BRST-symmetry leads to modified Slavnov-Taylor identities (mSTIs) for non-vanishing RG scales, k > 0 [2,8,[107][108][109][110][111][112]. The additional terms are generated by the BRST-variation of the regulator term and have a one loop form. They are similar in form and structure to the flow equation itself. The latter encodes the breaking or flow of scale invariance while the former encode the breaking or flow of BRST symmetry. The resulting mSTIs reduce to the standard STIs in the limit of vanishing RG scale, k → 0 , similar to the removal of the explicit breaking of scale invariance due to the regulator. Therefore, in both cases the underlying symmetry is restored in the limit of vanishing RG scale k → 0 . We emphasise that any regularisation scheme in momentum space leads to such a modification of BRST symmetry in terms of modified STIs. This is also well known from perturbation theory, where a cutoff regularisation, amongst other modifications, requires a gluon mass counter term in order to guarantee gauge invariance. Modified STIs are also present within other functional methods such as non-perturbative DSE and nPI approaches that rely on numerical momentum integrations.
To take the modification of the STIs of the vertices into account, we choose constant vertex dressings λc cA , λ A 3 and λ A 4 at the cutoff scale, k = Λ , such that the STIs for the running couplings,  (2) results [69,103].  are fulfilled at µ = 20 GeV, k = 0 . Here, the running couplings in (36) are obtained from the classical vertices, with the symmetric momentum configurationp . The mSTI of the gluon propagator implies a nonvanishing longitudinal gluon mass term at the cutoff scale [107]. In the perturbative regime, it can be shown that the transverse mass agrees with the longitudinal one, for details see [2]. However, while the longitudinal mass parameter vanishes at k = 0 , the transverse mass term encodes the gapping of the transverse gluon propagator at k = 0 . At the initial UV cutoff scale k = Λ the gluon mass parameter is uniquely determined by the mSTI and cannot be chosen freely. Its precise determination is at the root of confinement, which is encoded in the transverse mass gap at vanishing cutoff scale. Since the mass parameter is proportional to the cutoff, m 2 Λ ∝ α(Λ) Λ 2 , quadratic precision is required in its determination from the mSTI. The solution of this quadratic fine tuning problem requires both, a BRST-consistent quantitative level of the approximation, as well as sufficient numerical precision. Consequently, in truncated systems of flow equations, its computation from the mSTI at the required precision level is extremely challenging. The above brief discussion is detailed in [2]. Note that these statements hold also for other functional methods such as DSE and nPI approaches.
In the present work we utilise that it is possible to uniquely determine the gluon mass parameter by demanding a solution of the scaling type, for details see again [2]. We exploit that this also holds at finite temperature. Requiring scaling in the magnetic sector provides us with a unique value for the gluon mass parameter at each temperature. This procedure resolves the necessity of a BRST-consistent level of the approximation, but still requires quadratic precision in the fine-tuning. Further details are provided in App. G.

C. Numerical implementation
To solve the system of coupled flow equations, we use the tools established by the fQCD collaboration [23]. The tensorial flow equations are derived with DoFun [113]. Subsequently, the projected equations are traced with FormTracer [114], a Mathematica package that uses FORM [115][116][117][118] and has native support for finitetemperature applications. The output is exported as optimised C ++ code, which is then used within the computational framework of the fQCD collaboration. The latter uses the adaptive ordinary differential equation solver from the Boost library [119] and the adaptive multidimensional integration routine from [120], which implements [121,122].
In the derivation of the equations, tracing the fourgluon vertex equation, and in particular the gluon box diagrams, is the most challenging part. To this end, we generate FORM files with FormTracer [114] for each of the twelve permutations of the box diagrams. Executing one of these with FORM can take up to eight core days and intermediate expressions reach more than 1 TB in size. Since the resulting expressions are still very large, we sum all permutations, factorise all dressing functions and then use the simultaneous optimisation feature of FORM's optimisation routine [117] in combination with a parallelised version of FORM [123] to optimise the result. Concerning the numerical computation, integrating the flow once takes roughly one day on an ordinary quad-core desktop computer. This has to be done multiple times for each temperature due to the gluon mass parameter determination.

IV. RESULTS
The main results are displayed in Fig. 6 -11. We show results for the magnetic and electric dressing functions of propagators and vertices for various temperatures. For all correlators we find that the magnetic and electric dressings coincide for momenta p 2π T , and become degenerate with the vacuum dressings. This is required by the recovery of O(4) invariance. The convergence towards the vacuum dressings for small temperatures is explicitly checked in App. E. This apparently obvious property is actually non-trivial within frequency and momentum-dependent non-perturbative truncations.
We compare our gluon propagators to SU (2) [69, 103] and SU (3) [71] lattice results in Fig. 6 and 7. This comparison requires the setting of a relative scale as well as renormalisation, detailed in App. H. As a consequence, a potential relative offset of functional and lattice results has to be considered in addition to the systematic errors of the truncation, when juxtaposing the results from the different calculations. The comparison with SU (2) as well as SU (3) lattice data is legitimate because the truncation used in this work yields only a trivial dependence on the gauge group. This is the case because the colour traces can be taken without specifying the gauge group [114,124], and the only group constant appearing in the equations is the quadratic Casimir operator of the adjoint representation C A . Furthermore, C A occurs only in combination with the coupling at the renormalisation point α(µ) · C A ≡α(µ) . Thus, it can be absorbed into a redefinition of the running coupling, or, equivalently, the scale of the theory. Therefore, the propagators are identical for all groups, and the different couplings can be obtained by a trivial rescaling with C A . We emphasize that this is not a mere artefact of the approximation.
Perturbatively, the beta function of the pure gauge theory has a trivial group dependence up to three loops, see e.g. [125] for a recent discussion. See Sec. V for a discussion.
At low momenta, the electric and magnetic propagators show a qualitatively different behaviour. While the magnetic gluon propagator decreases almost monotonously with increasing temperature, the electric propagator increases at small temperatures. At high temperatures, where the growth of thermal contributions to the mass becomes dominant, also the electric gluon propagator decreases, see also Sec. II C and in particular Fig. 3. For the magnetic gluon propagator we find agreement with the lattice results on the 10 % accuracy level we expect from the truncation of the vertices. Furthermore, we see that the deviation takes its maximum for temperatures about the phase transition temperature, where we expect large-scale dynamical fluctuations to be most relevant. On the one hand, our truncation is tested maximally in this regime, and on the other hand discretisation and finite volume effects in the lattice calculation are strongest there. In contradistinction to the very satisfactory situation for the magnetic propagator, we observe a significant deviation about the phase transition temperature T c for the electric gluon propagator. However, the agreement is very good for small and large temperatures. As discussed in great detail in Sec. II B and V, the deviation about T c can be explained by the missing non-trivial A 0 -background in the present calculation.
The ghost propagator agrees qualitatively, but deviates quantitatively, from the lattice results, as shown in Fig. 8a. We discuss this point further in Sec. V. The ghost-gluon vertex is plotted in Fig. 8b. Interestingly, it is weaker around the phase transition temperature than in the vacuum. At high temperatures it shows a broader and less pronounced bump than at zero temperature.
The gluonic vertex dressing functions are shown in Fig. 9 and 10. The magnetic dressings of both vertices show scaling in the infrared. Contrarily, the correspond- ing electric components decouple at p ≈ 2π T and become constant in the infrared. We show the position of the zero crossing of the magnetic three-gluon vertex dressing function as a function of temperature in Fig. 11. At small temperatures the zero crossing moves towards lower momenta as the temperature is increased, since the three-gluon vertex is stronger at small and intermediate momenta for small temperatures, cf. Fig. 9. At high temperatures, the magnetic zero crossing rises linearly with the temperature. In contrast, the zero crossing of the electric three-gluon vertex dressing function disappears at T ≈ 40 MeV. Similarly, the electric dressing of the four-gluon vertex undergoes a drastic change from zero to small temperatures, where scaling is lost, and goes on to increase with growing temperature. At low momenta p 2π T the dimension of the theory is effectively reduced and the magnetic dressings behave as they do in three dimensions. In the case of the scaling solution, all magnetic dressing functions scale with a power-law [39,40,44], where 2n and m is the number of ghost and gluon legs, respectively. Due to dimensional reduction, the temperature-independent scaling coefficient κ is determined by three-dimensional Yang-Mills theory. Fitting the magnetic gluon propagators to (38) (3) . This agrees with the scaling exponent κ = 0.321(1) of the three-dimensional vacuum theory [100]. We summarise the different scaling coefficients in Tab. I.

V. DISCUSSION
In the previous section we have presented nonperturbative results obtained with the most comprehensive truncation within functional methods to date. The agreement of the magnetic propagator and the electric propagator for high temperatures is of the order of 10 % in the momentum regime relevant for hadronic observables. These small deviations can be attributed to lattice artefacts, the relative scale setting uncertainty, and the systematic error within our truncation. The latter stems from incomplete momentum dependencies of the vertices and missing non-classical tensors, see App. A and D for estimates of their respective importance. The electric propagator deviates from the lattice results at temperatures about the phase transition temperature. The explanation has already been indicated in Sec. IV and is discussed below.

A. Non-trivial backgrounds and their impact on electric and magnetic propagators
A potential source of the discrepancy of the electric gluon propagator near the phase transition temperature is an insufficient order in our approximation scheme. However, such deviations of the electric gluon propagator from lattice results were already observed in much simpler truncations [62]. Furthermore, if truncation artefacts were the main source, we would expect larger discrepancies also in the magnetic gluon propagator.
In contrast to this, the electric propagator, which is closely related to the order parameter L[ Ā 0 ], is partic- ularly sensitive to a non-vanishing background field [65]. As argued in Sec. II B, the non-trivial solution of the equation of motion,Ā 0 = 0 , is important in the temperature regime (25), that is This is exactly the temperature range where the deviations from the lattice results, which are evaluated on the equation of motion, are most pronounced. We expect a considerable improvement in the electric propagator if the correlation functions are evaluated on the non-trivial background. At this point, we want to emphasise that the observed deviations do not invalidate our results for the electric two-point correlator. It simply represents the correlation functions at a non-minimal configuration, cf. (18). Furthermore, these findings underline that Polyakov-enhanced low-energy effective models should be constructed inĀ 0 -backgrounds and the effective potential V [Ā 0 ] rather than Polyakov loop backgrounds and the Polyakov loop potential V [L] : the electric propagators agree on the 10 % level above T > ∼ 1.3 T c . This entails that the relevant background for the shifts in the Matsubara frequencies is Ā 0 .
The above analysis is also important for the discussion of the comparison of the present results with SU (2) and SU (3) lattice simulations. As discussed in detail in the last Sec. IV, the gauge group enters only at very high orders of the approximation in an expansion of the effective action around vanishing background. Thus, our results depend only trivially on the gauge group. However, the gauge group, and in particular the universality class, enters via the Polyakov loop background, or, more precisely Ā 0 . It has already been shown in [74,80] that the different orders of the phase transition for SU (2) and SU (N > 2) are encoded in the Polyakov loop potential V [Ā 0 ] and the respective expectation values ν in (24), rather than in the propagators. The Ising critical exponents for SU (2) are also extracted from critical fluctuations encoded in the effective potential, see [75] for Yang-Mills theory in Polyakov gauge and [102] for Landaugauge Yang-Mills theory. In [75,102] it also has been shown, that the critical fluctuations are the actual cause of the higher phase transition temperature in comparison to SU (N > 2) . Thus, the gauge group dependence of the order of the phase transition and the value of transition temperature are to leading order caused by the effective potential, and hence by the related expansion about the physical ground state, i.e. Ā 0 in the current setting.
We close the discussion of the propagators with the remark that the comparison of our results with the lattice results at small momenta p 2 Λ 2 QCD has to be taken with a grain of salt. The lattice results are of the decoupling type, while our results are of the scaling type. Consequently, possible non-perturbative gauge fixing effects have to be kept in mind, see e.g. [126][127][128][129]. This concerns in particular the ghost propagator, shown in Fig. 8a, which is more sensitive to the treatment of the Gribov copies than the gluon propagator [129].

B. Debye mass and perturbative regime
We find very good agreement of our non-perturbative Debye screening mass with two-loop hard thermal loop perturbation theory down to T ≈ 0.6 GeV, see Fig. 3. This remarkable agreement down to comparably low temperatures is in line with earlier findings, see e.g. [86][87][88][89][90][91][92][93][94]. In general, perturbative resummation schemes have been found to be applicable at surprisingly large couplings. An explanation of this unexpectedly large range of validity can be given by the structural similarity of higher order perturbative resummation schemes and the nonperturbative resummations performed within functional methods. This opens the door for applications of functionally assisted analytic perturbative computations beyond the validity bounds of perturbation theory, in particular to the transport and kinetic realm of heavy ion collisions.

C. Three-gluon vertex and its zero crossing
The magnetic three-gluon vertex dressing function has been studied on the lattice [72] and with a semiperturbative approximation of its DSE [73]. Both studies show a significant enhancement of the magnetic dressing at low momenta p ≈ 0.2 GeV for temperatures just below the critical temperature. While we also observe this effect qualitatively, see Fig. 9, we find a much weaker enhancement. This is consistent with the finding that our electric gluon propagator is weaker than the electric lattice propagator, cf. Fig. 7. This electric propagator enters the triangle diagram in the three-gluon vertex equation, which yields a positive contribution to the dressing function [53]. Thus, a stronger electric propagator increases the strength of the magnetic three-gluon vertex.
At zero temperature, the three-gluon vertex shows a zero crossing in four as well as in three dimensions [2, 6, 46, 49, 51-53, 59, 130-133]. Analytical studies show that it is caused by the divergent ghost triangle diagram. We find that the zero crossing persists in the magnetic dressing function for all temperatures. This stands in line with [73] but in contrast to [72], where the lowest investigated momenta show a positive sign at temperatures somewhat below the critical temperature. Here, we present an analytical argument for the persistence of the magnetic zero crossing at all temperatures. The argument is first presented for a vanishing gluonic background and is based on the infrared dominance of ghost loops. Finally we discuss the case of non-vanishing gluonic backgrounds relevant for temperatures about T c .
All gluonic diagrams are gapped below a certain scale, whereas the ghost triangle effectively behaves like the corresponding three-dimensional diagram for p 2 (2π T ) 2 . Therefore, it causes a divergence in the magnetic threegluon vertex dressing function at low momenta for all temperatures, and thus, the magnetic zero crossing cannot vanish. At high temperatures, this zero crossing moves then to higher scales, which is in line with the high temperature limit and [72]. This qualitative argument is actually independent of the type of the solution, since the three-dimensional ghost triangle diagram diverges with a power-law in the case of the scaling solution and linearly [6,49,51] in the case of the decoupling solution. We find that the zero crossing of the electric component vanishes at a temperature of T ≈ 40 MeV. This can be understood by observing that the zero mode of the ghost triangle diagram, evaluated at zero external Matsubara frequencies, contributes to the magnetic three-gluon vertex dressing, but vanishes analytically if projected with the electric three-gluon vertex projection operator. Our numerical results show precisely the expected behaviour, see Fig. 11 and 9.
We extend the argument to the case of non-vanishing backgrounds. They introduce a colour structure in the ghost propagator and the ghost-gluon vertex. After diagonalisation, we are left with gapped and ungapped modes in the ghost propagator, as well as background- dependent and background-independent (colour) tensor structures in the ghost-gluon vertex. The remaining ungapped ghost modes couple to the latter tensor structure, which is nothing but the original tensor structure at vanishing background. Therefore, the background simply leads to a weakening of the infrared dominance by gapping some, but not all, ghost modes. Accordingly, the zero crossing moves towards smaller momenta, but does not disappear, in the presence of non-trivial backgrounds.
Furthermore, for small temperatures T /Λ QCD → 0 , the gapping of the ghost occurs only at very small momenta p 2 < ∼ (2π T ) 2 , and we are left with the temperature regime (25), in which a weakening of the infrared ghost dominance is to be expected. This structure is compatible with the results in [72], where no zero crossing was observed at temperatures about T c in the accessible momentum regime. In our opinion, it would therefore be interesting to extend the analysis of [72] to smaller momenta.

VI. SUMMARY & OUTLOOK
We have presented non-perturbative first-principles results for the finite-temperature Landau-gauge Yang-Mills correlation functions, obtained from the functional renormalisation group. Our comprehensive truncation of the effective action includes the computationally especially expensive magnetic and electric components of the purely gluonic vertices. We gauged our truncation by comparing to propagator results obtained in lattice simulations and found very good agreement for the magnetic gluon propagator. Our result for the Debye screening mass shows excellent agreement with two-loop hard thermal loop perturbation theory at high temperatures and the electric gluon propagator compares very well to lattice results for all temperatures except T ∈ (0.5 T c , 1.3 T c ) . We have argued that the deviations in this regime are related to the different backgrounds used. Particular focus was also put on the fate of the zero crossing in the three-gluon vertex at finite temperature. In the electric component of the three-gluon vertex we found the disappearance of the zero crossing at a very small temperature. The magnetic zero crossing also moves towards lower momenta for small temperatures but it never vanishes. At high temperatures, its position increases linearly with the temperature. We gave an analytic argument for the observed qualitative behaviour of the zero crossing in the magnetic and electric components.
The presented first-principles results for the finitetemperature correlation functions of Yang-Mills theory form the foundation for a number of subsequent studies. First and foremost, the capability to perform nonperturbative first-principles studies of gauge theories at finite temperatures provides a crucial prerequisite for the investigation of the QCD phase structure. In particular, combining the advancements of this work with those of a recent calculation of the correlators of two-flavour QCD [3], will allow us to investigate the properties of quantum chromodynamics at finite temperature and density from first principles. Furthermore, the presented correlators can be used to compute thermodynamic quantities like the pressure, the shear viscosity, as well as the Polyakov loop potential and spectral functions, the latter being notoriously difficult to obtain. Additionally, it is suggestive to use the remarkable agreement of fully nonperturbative results with resummed perturbative results, in particular in the Debye mass, to devise functionally assisted analytic applications for the transport and kinetic regime in heavy ion collisions. Finally, we expect that improving the current investigation by including a nonvanishing background field and non-vanishing Matsubara modes will lead to the disappearance of the discrepancy in the electric gluon propagator near the phase transition temperature.

ACKNOWLEDGMENTS
We thank Lukas Corell, Ouraman Hajizadeh, Markus Q. Huber, Axel Maas, Richard Williams, and Nicolas Wink for discussions. This work is supported by EMMI, the grant ERC-AdG-290623, the FWF through Erwin-Schrödinger-Stipendium No. J3507-N27, the Studienstiftung des deutschen Volkes, the DFG through grant STR 1462/1-1, the BMBF grant 05P12VHCTG, and is part of and supported by the DFG Collaborative Research Centre "SFB 1225 (ISOQUANT)". It is also supported in part by the Office of Nuclear Physics in the US Department of Energy's Office of Science under Contract No. DE-AC02-05CH11231.  [2] and SU (3) lattice data [134]. The lattice results are renormalised as in [2]. Newer lattice results [135] agree with [134] if the largest physical volumes are compared. .

Appendix A: Regulator and truncation dependence
The regulators in (31) are parametrised by where we dress the regulators withZ M A,k andZ c,k as in [2]. Since Π ⊥ µν (p) = Π M µν (p) + Π E µν (p) , (A1) implies the same regularisation for electric and magnetic modes. Due to its advantages for the evaluation of thermodynamic quantities [62,97] we use the exponential regulator shape function, with m = 2 . This is in contrast to the vacuum calculations in [2], which were performed with a smoothed version of the flat regulator [105].
In Fig. 12, we show vacuum results obtained with the flat and the exponential regulator. Clearly, the results obtained with the symmetric-point approximation, defined by (13), agree very well. However, they show a higher bump than the lattice results. This is due to the symmetric-point approximation used for the vertices. This discrepancy vanishes if more momentum dependencies are included as shown in Fig. 12, cf. [2] for a thorough discussion. An extension of the current finitetemperature investigations beyond the symmetric-point approximation is deferred to future work.

Appendix B: Tensor Splitting
At vanishing Matsubara frequencies, not all tensors that are obtained by contracting the classical tensor structures with all possible combinations of electric and magnetic projectors are linearly independent, see App. C. Since we calculate the dressing functions at vanishing Matsubara mode, we can compute only a restricted set of dressings. This implies that we have to approximate the remaining degenerate dressing functions from this reduced set of dressings in order to obtain the correct UV behaviour and to recover the vacuum results in the zerotemperature limit. We use for the ghost-gluon vertex, for the three-gluon vertex and [Γ A 4 ] abcd µνρσ (p, q, r) = [S (4) for the four-gluon vertex. Here and in the following we leave the momentum arguments implicit. Although the dressing functions of some tensors coincide in our approximation, we explicitly show the splitting to make the construction of the approximation apparent. The magnetic dressing functions appear in more than one tensor structure, and we evaluate them by projecting onto the purely magnetic tensors structure for every vertex, see App. D.
Due to the O(4)-symmetry of the vacuum, this approximation becomes exact for large momenta p 2 (2π T ) 2 , which are not affected by finite-temperature effects. In the limit of vanishing Matsubara frequencies, the dimension of the tensor space is reduced, see App. C. Therefore, and for the four-gluon vertex, Thus, for p 2 (2π T ) 2 the classical vertex dressings are fully described by the remaining basis tensors, to wit, The tensor bases for the propagators as well as for the ghost-gluon vertex are complete, and therefore the projection onto the dressings is unique. For the gluonic vertices we do not take the full transverse tensor bases into account. Consequently, already in the vacuum, any projection is an approximation that relies on the assumption that non-included basis elements are small. If the flows are projected onto their electric and magnetic components, the incompleteness of the bases can lead to intricate complications. The reason is that the magnetic and electric projection operators can yield differing contributions from non-classical tensor structures that are created by quantum fluctuations. As an immediate consequence, the magnetic and electric dressings differ then by momentum dependent terms. This effect occurs already at vanishing temperature, and is therefore in contradiction with the O(4)-symmetry of the vacuum. If one uses a complete basis, projecting with magnetic and electric projection operators does not spoil the O(4)-symmetry although the projection operators themselves are not O(4)-symmetric.
In the following two subsections we discuss in detail the quantitative relevance of these effects caused by the incomplete bases for the gluonic vertices. In order to disentangle genuine finite-temperature contributions from these projection artefacts, we consider only vacuum flows. By splitting the projection into electric and magnetic components and comparing them to the O(4)symmetric projection, we are able to quantify these basis artefacts. Unfortunately, we find that the emergence of certain non-classical tensors yields sizeable artefacts on the dressing of the classical tensor structure of the four-gluon vertex. As discussed in detail in this and the following two App. E and F, implementing a proper treatment of these artefacts of the incomplete bases turns out to be vital to obtain the correct UV behaviour and cutoff independence of the finite-temperature results.

Three-gluon vertex
We project onto the magnetic and electric components of the three-gluon vertex by as generalisation of the vacuum projection In explicit numerical checks we find that the projections (D1) and (D2) agree to the per mille level at T = 0 and therefore also for k 2π T . We conclude that our projection is not sensitive to the possible emergence of non-classical tensors structures in the three-gluon vertex. This is also consistent with the sub-leading importance of non-classical tensor structures found in earlier threegluon vertex studies [53].

Four-gluon vertex
We project onto the vacuum dressing function with Assuming vanishing non-classical tensor structures, this generalises to for the magnetic and the electric components. If the only tensor generated by the flow were the classical one, the projections (D3) and (D4) would yield λ A 4 = λ M A 4 = λ E A 4 . However, this equality can be spoiled by the presence of non-classical tensors, which are in general created by the flow equation. Consider, for example, the following O(4)-and Bose-symmetric non-classical tensor: [Γ (4) A 4 ,ncl ] abcd µνρσ (p, q, r, s) = f abn f cdn · (q + s) µ (q + s) ρ (p + r) ν (p + r) σ + perm. . (D6) Inserting (D6) into (D3) and (D4) yields differing contributions to the dressing functions λ A 4 , λ M A 4 , and λ E A 4 . Therefore, O(4)-invariance is lost due to the incompleteness of the basis that was used to construct the projection operators, (D3) and (D4).
In Fig. 13 we show the vacuum flows of the four-gluon vertex obtained with different projection operators and identical vacuum vertices on the right hand side of the flow equation. In contrast to the three-gluon vertex, we find a considerable difference in the resulting momentum dependence of the projections (D3) and (D4). We conclude, that sizeable non-classical tensors, which affect the difference between the magnetic and electric projection operators, are generated. As an immediate consequence, the O(4)-symmetric limit at T → 0 is spoiled by the presence of these tensors since λ M A 4 (T = 0) = λ E A 4 (T = 0) . A simple estimate of the unphysical projection artefacts of these non-classical tensors is given by the vacuum differences of the projections (D3) and (D4), which are also shown in Fig. 13. Assuming that this unphysical difference depends only mildly on the temperature, a natural strategy to account for this artefact is to subtract (D7) from the finite-temperature flows. However, there is an additional complication in the case of the scaling solution. The vertex dressings obey a power law behaviour at small momenta, see (38), and the corresponding exponent changes as one goes from the vacuum to finite temperature. Although this behaviour of the correlators at very small momenta does not affect any observables, it has to be taken into account when subtracting (D7) from the finite-temperature flows. Conse- quently, we modify the flows of the magnetic and electric components by The purpose of the smoothed step function, θ (k, k c ) = 1 is to provide a transition from the corrected flows to the pure finite-temperature flows with the correct scaling behaviour at very low momenta. We set the transition scale k c to k c = min (λ c 2π T, Λ c ) , which is defined in terms of the parameters Λ c and λ c . The modified dressings fulfill This guarantees that we recover the vacuum results in the limit of vanishing temperature while our best estimates for the basis artefacts are subtracted above the transition scale k c . In order to investigate the influence of the precise value of the transition scale, we vary Λ c and λ c in reasonable ranges. Since temperature effects are expected to be small at momentum scales k ≥ 2π T , λ c should be of order unity and we vary it from 1 to 2 . Furthermore, the gapping scale of the gluon propagator gives us an estimate on the scale below which no phenomenologically important effects are to be expected. Consequently, we vary Λ c between the location of the maxima of the gluon propagator and the gluon propagator dressing, i.e. Λ c ∈ [0.3, 1] GeV. We find only a mild (10 %) dependence of the four-gluon vertex dressings on these parameters as shown in Fig. 14. Since the four-gluon vertex is the least important of all classical tensors in the self-consistently coupled system, we find that the dependence of all other dressings on these parameters is even smaller. For example, the induced uncertainty on the electric gluon propagator is at most 3 %, but for a wide range of temperatures and momenta it is even smaller than 0.5 %. In all cases the dependence on the smoothing parameter, = 0.05 , is negligible.  This property enables us to reduce the computational effort by one to two orders of magnitude. First we compute the T = 0 effective action, starting at a large perturbative scale of typically k = Λ = 60 GeV from the classical action. To obtain the temperature-dependent 1PI effective action, we integrate the flow equation starting from the effective average action Γ Λ T (T ) = Γ Λ T (T = 0) at a lower, temperature-dependent cutoff scale, Here, Λ min T has been introduced to avoid the interference of the lowered starting scale with the dynamical mass generation of the gluon. This is necessary, because the scaling condition forces us to modify the input at Λ T by a gluon mass term, see App. G. Consequently, we choose Λ min T as the scale where the vacuum gluon propagator dressing becomes maximal, i.e. Λ min T = 0.955 GeV. We show the dependence of the longitudinal gluon propagator on the physical start scale λ in Fig. 15b. The gluon propagator as well as all other quantities do not depend on the start scale for λ ≥ 2 . In our numerical computation we use λ = 4 although λ = 2 is sufficient, as argued in App. G.
To demonstrate the advantage of the temperaturedependent initial scale, we consider the numerical vacuum integration, Numerically, it is advantageous to choose a k-dependent numerical cutoff L = l k , where l = 3 is sufficient for the exponential regulator due to the regulator derivative appearing in all diagrams. This persists in the Matsubara formalism and we can limit the summation to frequencies ω = 2π T n smaller than L = lk , (F4) Thus, the number of required integrand evaluations grows linearly with k as well as with 1/T . For small temperatures, the increased number of evaluations due to the growing number of small Matsubara modes is therefore compensated by the shrinking initial scale, at least down to T = Λ min T /(2πλ) .   (26). The dashed lines are fits of (I1) with c a = 0 , i.e. (27), to the large points. Small points show G E T (x) beyond the fit regions. The fitted screening masses as a function of temperature are shown in Fig. 3. .
The SU (3) lattice results from [71] do not include the vacuum case T = 0 . Therefore, we allow for a scale mismatch by introducing a temperature-independent relative scale factor r s , in addition to the temperatureindependent wave function renormalisation constant z L . We determine r s and z L by fitting the magnetic gluon dressing function, 1/ z L Z M A (r s T , r s p) , simultaneously for all temperatures to all lattice points above p ≥ 0.5 GeV. Subsequently, we use r s and z L to rescale the magnetic as well as the electric lattice propagators to our data. We find the relative scale mismatch r s − 1 to be small, of the order of 2 %. The temperatures in [71] are given in units of GeV. In order to simplify the discussion, we convert the temperatures into units of the critical temperature, using their value for the SU (3) phase transition temperature, T c = 270 MeV.